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Introduction 


A  novel  imaging  technology,  scanning  microwave-induced-acoustic  tomography,  will  be 
developed  for  breast  imaging.  X-ray  mammography  and  ultrasonography  are  the  current  clinical 
tools  for  breast-cancer  screening  and  detection.  Mammography  is  the  “gold  standard”,  however, 
uses  ionizing  radiation  and  has  difficulties  imaging  pre-menopausal  breasts,  which  are 
radiographically  dense.  Ultrasonography  is  an  adjunct  tool  to  x-ray  mammography  and  cannot 
detect  many  of  the  nonpalpable  tumors.  The  cure  rate  of  breast  cancers  is  improved  if  they  are 
detected  early.  To  provide  a  new  non-invasive,  non-ionizing  diagnostic  tool  for  detection  of  early 
breast  cancers,  we  will  develop  real-time  microwave-induced-acoustic  tomography  for  breast 
imaging.  Microwave-induced-acoustic  tomography  is  based  on  the  photoacoustic  effect, 
generation  of  acoustic  wave  by  deposition  of  short-pulse  electromagnetic  energy  safely  into 
biological  tissues.  The  microwave  for  this  technology  is  short-pulsed,  and  its  power  is  within  the 
IEEE  safety  limits.  The  microwave-induced  acoustic  wave  is  then  detected  with  an  ultrasonic 
detector  for  imaging.  The  contrast  between  tumors  and  normal  tissues  in  the  microwave  regime  is 
significantly  better  than  other  imaging  modalities.  Cancerous  breast  tissues  are  found  to  be  2-5 
times  more  strongly  absorbing  than  surrounding  normal  breast  tissues  in  the  microwave,  which 
has  been  attributed  to  an  increase  in  bound  water  and  sodium  within  malignant  cells.  However, 
pure-microwave  imaging  is  fundamentally  limited  to  poor  resolution  (on  the  order  of  10  mm) 
because  of  the  large  wavelength  of  microwave.  Ultrasonic  imaging  has  good  resolution  (on  the 
order  of  1  mm)  but  has  a  poor  contrast  between  tumors  and  normal  tissues.  Microwave-induced- 
acoustic  tomography  combines  the  contrast  advantage  of  pure-microwave  imaging  and  the 
resolution  advantage  of  pure-ultrasonic  imaging,  therefore,  has  the  potential  for  detection  of  early 
breast  cancers  and  for  assessing  and  monitoring  treatments  as  well. 


Body 

In  this  section,  we  present  our  study  of  pulsed-microwave-induced  thermoacoustic  tomography  in 
biological  tissues.  A  short-pulsed  microwave  source  was  used  to  irradiate  the  tissue  samples,  and 
the  thermoacoustic  waves  excited  by  thermoelastic  expansion  were  then  measured  by  a  wide¬ 
band  ultrasonic  transducer  along  a  circular  path  that  encloses  the  sample  under  study.  The 
acquired  data  were  then  used  to  reconstruct  the  microwave  absorption  distribution.  Both  an 
exact  reconstruction  solution  and  an  approximate  modified  backprojection  algorithm  were 
derived.  Experiments  demonstrated  that  the  images  calculated  by  the  backprojection  method 
agreed  with  the  original  samples  very  well,  and  the  spatial  resolution  in  reconstruction  was  as 
good  as  0.5  millimeter  (500  micrometers). 

Introduction  to  thermoacoustic  tomography 

In  thermoacoustic  tomography,  a  short-pulsed  microwave  source  is  used  to  irradiate  the  tissue. 
Absorbed  microwave  energy  causes  thermoelastic  expansion  and  radiates  thermoacoustic  waves 
from  within  the  irradiate  tissue.  The  relatively  long  wavelength  of  the  microwave,  e.g.,  ~3  cm  at 
3  GHz  in  tissues,  serves  to  illuminate  the  tissue  homogeneously.  The  microwave  heating  must 
be  rapid  to  produce  thermoacoustics  waves;  in  other  words,  static  temperature  distribution  or 
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slow  heating  cannot  produce  thermoacoustic  waves.  A  wide-band  ultrasonic  transducer  can  then 
be  employed  to  acquire  the  thermoacoustic  signals  excited  by  thermoelastic  expansion,  which 
carries  the  microwave  absorption  property  of  the  tissue.  The  ultrasonic  transducer  is  very 
sensitive  in  detecting  small  vibrations  from  an  object  that  are  caused  by  weak  energy  absorption. 

The  key  problem  with  this  technique  is  how  to  determine  the  microwave  absorption 
distribution  from  the  measured  data,  i.e.,  how  to  map  the  inhomogeneity  of  the  tissue.  One 
approach  is  to  use  focused  ultrasonic  transducers  to  localize  the  thermoacoustic  sources  in  linear 
or  sector  scans  and  then  construct  the  images  directly  from  the  data  as  is  often  done  in  pulse-echo 
ultrasonography.  An  alternative  method  is  to  use  wide-band  unidirectional  point  detectors  to 
acquire  thermoacoustic  data  and  then  reconstruct  the  microwave  absorption  distribution.  To 
date,  we  have  not  seen  an  exact  inverse  solution  for  this  specific  problem,  although  some 
researchers  have  arrived  at  approximate  reconstruction  algorithms,  such  as  the  weighted  delay- 
and-sum  method,  the  optimal  statistical  approach,  and  the  Radon  transform  in  far  field 
approximation. 

Based  on  spherical  harmonic  functions,  we  first  deduced  an  exact  solution  of  the  problem 
in  the  three-dimensional  case,  which  can  be  carried  out  in  the  frequency  domain.  We  assume  that 
the  wide-band  unidirectional  ultrasonic  transducer  is  set  on  a  spherical  surface,  which  encloses 
the  sample  under  investigation.  The  data  acquired  from  different  directions  are  sufficient  to 
allow  us  to  reconstruct  the  microwave  absorption  distribution.  In  our  case,  the  diameter  of  the 
sphere  of  detection  is  much  larger  than  the  ultrasonic  wavelength.  Next,  an  approximate 
algorithm  is  deduced,  which  is  a  modified  backprojection  of  a  quantity  related  with  the 
thermoacoustic  pressure.  This  approximate  algorithm  can  be  carried  out  in  the  time  domain  and 
is  much  faster  than  the  exact  solution.  We  have  also  tested  a  set  of  tissue  samples.  These 
experiments  demonstrate  that  the  images  calculated  by  the  modified  backprojection  method  agree 
with  the  original  samples  very  well.  Moreover,  the  images  have  both  the  high  contrast  associated 
with  pure-microwave  imaging  and  the  0.5-millimeter  spatial  resolution  associated  with  pure- 
ultrasound  imaging. 

Theory  of  thermoacoustic  tomography 

Fundamentals  ofthermoacoustics 

Thermoacoustic  theory  has  been  discussed  in  many  literature  reviews  such  as.  Here,  we  briefly 
review  only  the  fundamental  equations.  If  the  microwave  pumping  pulse  duration  is  much 
shorter  than  the  thermal  diffusion  time,  thermal  diffusion  can  be  neglected;  consequently,  the 
thermal  equation  becomes 

pCp^-T(r,t)  =  H(r,t),  (1) 

ot 

where  pis  the  density,  Cp  is  the  specific  heat,  T(r,t)  is  the  temperature  rise  due  to  the  energy 
pumping  pulse,  and  H(r,t)  is  the  heating  function  defined  as  the  thermal  energy  per  time  and 
volume  deposited  by  the  energy  source.  We  are  interested  in  tissue  with  inhomogeneous 
microwave  absorption  but  a  homogeneous  acoustic  property.  The  two  basic  acoustic  generation 
equations  in  a  homogeneous  medium  are  the  linear  inviscid  force  equation 

u(r,  0  =  ~vP(r>  0  (2) 
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and  the  expansion  equation 


V  •  u(r,  t)  =  +  pT  (r,  0  (3) 

P° 

where  /?  is  the  isobaric  volume  expansion  coefficient,  cis  the  sound  speed,  u(r,/)  is  the 
acoustic  displacement  and  p(r,t )  is  the  acoustic  pressure. 

Combining  the  above  three  equations,  the  pressure  p(r,t)  produced  by  the  heat  source 
H(r,t )  obeys  the  following  equation 


(4) 

The  above  equation  is  a  typical  scalar  Helmholtz  equation.  The  solutions  based  on  Green's 
function  can  be  found  in  the  literature  of  physics  or  mathematics.  A  general  form  can  be 
expressed  as 


P(r>t)  = 


£ 

4  nCp 


rrr  d3r  dH(r',t') 
—  r'  dt' 


(5) 


The  heating  function  can  be  written  as  the  product  of  a  spatial  absorption  function  and  a  temporal 
illumination  function: 


H(r,t)  =  A(r)I(t). 


(6) 


Thus,  p(r,t )  can  be  expressed  as 

P(r,t) 


£ 

4 *Cp 


A(r')l'(t') . 


(7) 


Exact  re construction  theory 

We  first  solve  the  problem  where  the  pulse  pumping  is  a  Dirac  delta  function  as 

I{t)  =  I0S(t).  (8) 

Suppose  the  detection  point  on  the  spherical  surface  r  =  r0 ,  which  encloses  the  sample  (Fig.  1). 
By  dropping  the  primes,  the  pressure  equation  may  be  written  as 

P(.r0,t)  =  V \\\^rA{r)  ^  ’  (9) 

where  n  =  .  The  inverse  problem  is  to  reconstruct  the  absorption  distribution  A( r)  from  a 

C, 

set  of  data  p(rQ,t)  measured  at  position  r0 .  Taking  the  Fourier  transform  on  variable  t  of  the 


above  equation  and  denoting  k  =  — ,  we  get 

c 


p(r0,a))  =  -icon  fJ^3^(r)eX^r°_r|Ii)- , 


(10) 


where  following  Fourier  transform  pair  exists: 
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p(r0,a))  =  j/?(r0,0  exp(iax)dt , 
p(r0,t)  =  —  \p(r0,co)exp(-icot)dG). 


(11) 


(12) 


The  exact  inverse  solution  to  the  pressure  equation  (Eq.  9)  can  be  derived  on  the  basis  of  the 
spherical  harmonic  function: 


(2m  +  l)jjkr) 


^■(n  •“<>)’ 


(13) 


where  n  =  r/r,  n0=r0/r0,  j,(-)  and 

/?/<l)(-)  are  the  spherical  Bessel  and 

Hankel  functions,  respectively;  P{) 
represents  Legendre  polynomial.  This 
inverse  solution  involves  summation 
of  a  series  and  may  take  much  time  to 
compute.  Therefore,  it  is  desirable  to 
further  simplify  the  solution. 

Modified  backoroiection 
In  the  experiments,  the  detection 
radius  r0  is  much  larger  than  the  wavelengths  of  the  thermoacoustic  waves  that  are  useful  for 

imaging.  Because  the  low-frequency  components  of  the  thermoacoustic  signal  do  not 
significantly  contribute  to  the  spatial  resolution,  they  can  be  removed  by  a  filter.  Therefore,  we 
can  assume  |&|r0  » 1  and  use  the  asymptotic  form  of  the  Hankel  function  to  simplify  the  above 
exact  inverse  solution.  The  following  two  identities  are  involved: 
exp(-/£|r0  -r|)  _-ik 


Fig.  1.  Acoustic 
detection  scheme. 
The  ultrasonic 
transducer  at 
position  r0  records 
the  thermoacoustic 
signals  on  a 
spherical  surface 
with  radius  |r-rn|. 


When  |&|r0  » 1 , 


4/r|r0-r| 
hmfa o)" 


£  (2m  +  \)jm  (kr)h(m2)  (kr0  )Pm  (n  •  n0 ) ; 

4*  to 


1 


h^(kr0) 


1 


(^o)2 


■+d 


i 


.(^o)4 


(14) 


(15) 


JJ 


where  hj2)(-)  is  the  spherical  Hankel  function  of  the  second  kind.  After  some  mathematical 
operations,  the  approximate  inverse  solution  becomes 


A(r)  =  - 


2rn]C  £ 


dt 


(16) 


i.e.. 
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The  above  equation  shows  that  the  absorption  distribution  can  be  calculated  by  means  of 

backprojection  of  the  quantity  instead  of  the  acoustic  pressure  itself.  This 

v  t  dt  ,  M 

c 

approximate  algorithm  involves  less  computing  time  than  the  exact  inverse  solution. 

For  initial  investigations,  we  reconstruct  individual  cross  sections  of  samples.  In  these 
cases,  the  backprojection  is  carried  out  in  a  circle  around  the  cross  sections,  and  the  approximate 
inverse  solution  can  be  simplified  as 

m  =  — Lr  (18) 

June  •  t  dt  ,_k±l 


Experimental  method  of  thermoacoustic  tomography 

Experimental  setu]) 

Fig.  2  shows  the  experimental  setup.  A  plexiglass  container  is  filled  with  mineral  oil.  A  rotation 
stage  and  an  unfocused  ultrasonic  transducer  are  immersed  inside  it  in  the  same  x-y  plane.  The 
slice  sample  can  be  put  in  the  rotation  stage  horizontally.  The  transducer  points  horizontally  to 
the  rotation  center  and  detects  the  acoustic  signal  from  the  sample.  A  step  motor  directly  drives 
the  rotation  stage  while  the  transducer  is  fixed.  Obviously,  this  is  equivalent  to  having  a 
transducer  rotationally  scanning  the  sample.  The  transducer  (V323,  Panametrics)  has  a  central 

frequency  of  2.25  MHz  and  a  diameter  of  6  mm.  _ _ _ — ===n 

The  microwave  pulses  are  step  motor  - _ _ _  Driver - • 

transmitted  from  a  3 -GHz  n  Ur-1  cjp  n  computer 

microwave  generator  with  a  energy  "  ~  '  ir  — *  Jj 

of  10  mJ/pulse  and  a  pulse  width  of  ■=□ - L>=====0:=  ~  I - 

0.5  |is.  A  function  generator  I  1  I - J 

(Protek,  B-180)  is  used  to  trigger  Coupling .  Pulse  amplifier 

the  microwave  generator,  control  4  II 

its  pulse  repetition  frequency,  and  &-►  u,  u.J  trp-l  - 

synchronize  the  oscilloscope  y  sample  Transducer  Function  Oscilloscope 

J  _  .  ii  — i  qenerator 

sampling.  Microwave  energy  is  n  T 

delivered  to  the  sample  from  below  1 - 

by  a  rectangular  waveguide  with  a  U| — I - 1 — ,  U  II  I - 1  I - 1 

cross  section  of  72  mm  X  34  mm.  Microwave  generator  >j -  I - 1 

A  personal  computer  is  Fig.  2.  Experimental  setup, 

used  to  control  the  step  motor  in 

rotating  the  sample.  The  signal  from  the  transducer  is  first  amplified  through  a  pulse  amplifier, 
then  recorded  and  averaged  200  times  by  an  oscilloscope  (TDS640A,  Tektronix),  and  finally 
transferred  to  a  personal  computer  for  imaging. 

Lastly,  we  want  to  point  out  that,  in  our  experiments,  the  distance  r0  between  the  rotation 

center  and  the  surface  of  the  transducer  is  4.3  cm.  In  the  frequency  domain  (60  KHz-1.8  MHz), 
\k\r0=2 nrjlc  with  1.5  mm/ps,  we  get  10  <  |/c|r0  <  330 .  Therefore,  the  required  condition 

\k\r0  » 1  for  the  modified  backprojection  algorithm  is  satisfied. 


"Couplfng 

medium 


Sample  Transducer 


Microwave  generator « 


Pulse  amplifier 


Function 

generator 


Oscilloscope 


Fig.  2.  Experimental  setup. 
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Technical  consideration 

The  ultrasonic  transducer  is  not  a  real  point  detector.  For  simplicity,  we  can  ignore  its  size  if  we 
put  it  far  away  from  the  sample.  However,  we  still  have  to  consider  the  impulse  response  R(t) 
of  the  transducer  and  the  pumping  duration  /(/)  of  the  microwave  pulse.  In  general,  the 
measured  piezoelectric  signal  can  be  written  as  a  convolution: 

S(r0,t)  =  p(r0,t)*I(t)*R(t),  (19) 

where  p(r0,t)  is  the  thermoacoustic  signal  with  delta-pulse  microwave  pumping.  In  the 

frequency  domain,  the  above  equation  can  be  written  as 

S(r0,©)  =  /*r0, ©)/(©)*(©),  (20) 

where 

/(to)  =  J/(/)exp(/o>/)dr ,  (21) 

— OO 

+oo 

R((£>)  =  ^R(t)exp(iQ)t)dt .  (22) 


Therefore,  d/fotnO  can  be  calculated  by  an  inverse  Fourier  transformation, 
dt 


MTqJ)  _ 
dt 


FFT~'<!  l0]S(r^  F{0))\ 


I(0))R(co) 
1  +7-icoS(r0,O)) 


(23) 


-if 


I{o))R((o) 


F(cd)  exp(-icot)dt 


where  F(a>)  is  a  wide  band-pass  filter,  which  is  used  to  eliminate  the  noise  at  high  frequencies 
as  well  as  the  low  frequency  component  to  guarantee  the  condition  \k\rQ  » 1  for  the  modified 
backprojection. 

In  our  experiments,  I(t)  is  approximately  a  rectangular  function  with  duration  T  =  0.5  ps . 
The  ultrasonic  transducer  is  of  the  videoscan  type  with  a  central  frequency  fo  =  2.25  MHz.  The 
generated  thermoacoustic  signal  mainly  exists  in  a  frequency  range  below  1 .8  MHz.  Therefore,  a 
band-pass  filter  F(o))  may  be  employed  in  data  processing,  which  lets  only  the  signal  in  the 
range  between  60  KHz  and  1 .8  MHz  pass  through. 


Results  and  discussion  of  thermoacoustic  tomography 

Image  contrast 

Image  contrast  is  an  important  index  for  biological  imaging.  Fig.  3(a)  shows  a  tested  sample, 
which  was  photographed  after  the  experiment.  The  sample  was  made  according  to  the  following 
procedure.  First,  we  cut  a  thin  piece  of  homogeneous  pork  fat  tissue  and  shaped  it  arbitrarily  to 
form  a  base.  Its  thickness  is  5  mm  and  its  maximum  diameter  is  4  cm.  Then  we  used  different 
screwdrivers  to  carefully  make  two  pairs  of  holes  that  were  approximately  4  mm  and  6  mm  in 
diameter,  respectively.  Finally,  one  big  and  one  small  hole  on  the  left  side  was  filled  with  pork 
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Fig.  3.  (a)  Photograph  of 

the  cross-section  of  a  tissue 
sample;  (b)  Reconstructed 
image. 


muscle,  while  the  big  and  small  hole  on  the  right  side  were  filled  with  pork  fat  of  the  same  type 
as  that  which  made  up  the  base^ 

In  the  experiment, 
the  transducer  rotationally 
scanned  the  sample  from  0 
to  360  degrees  with  a  step 
size  of  2.25  degrees.  We 
used  the  160  series  of  data 
to  calculate  the  image  by 
our  modified 

backprojection  method. 

The  reconstructed 
image  is  shown  in  Fig.  3(b). 

The  outline  and  size  of  the 
fat  base  as  well  as  the  sizes 
and  locations  of  the  two 
muscle  pieces  are  in  good 
agreement  with  the  original  sample  in  Fig.  3(a).  The  high  contrast  is  due  to  the  low  microwave 
absorption  capacity  of  fat  and  the  high  absorption  capacity  of  muscle:  at  3  GHz,  the  penetration 
depth  for  muscle  and  fat  are  1 .2  cm  and  9  cm,  respectively.  The  two  pieces  of  fat  are  not  visible 
in  the  image  Fig.  3(b),  which  means  the  minute  mechanical  discontinuity  between  the  boundaries 
of  muscle  and  fat  does  not  contribute  much  to  the  thermoacoustic  signal.  On  the  contrary,  the 
discontinuity  improves  the  strength  of  the  echo  sounds  in  pure-ultrasound  imaging. 


Horizontal  axis  x  (mm) 


(b) 


-0.6  -0.4  -0.2  -0.0 


SjmtUU  resolution 

Spatial  resolution  is  another  important  index  for  biological  imaging.  We  used  samples  with  a  set 
of  small  thermoacoustic  sources  to  test  the  resolution.  One  tested  sample  is  shown  in  Fig.  4(a), 
which  was  also  photographed  after  the  experiment  was  completed. 

The  sample  was  made  according  to  the  following  procedure.  First,  we  cut  a  thin  piece  of 
homogeneous  pork  fat  tissue  and  made  it  into  an  arbitrary  shape.  Its  thickness  was  5  mm  with  a 
maximum  diameter  of  4  cm.  Then  we  used  a  small  screwdriver  to  carefully  make  a  set  of  small 
holes  about  2  mm  in  diameter.  In  the  meantime,  we  prepared  a  hot  solution  with  5%  gelatin, 
0.8%  salt  and  a  drop  of  dark  ink  (to  improve  the  photographic  properties  of  the  sample).  Next, 
we  used  an  injector  to  inject  a  drop  of  the  gelatin  solution  into  each  small  hole  and  subsequently 
blew  out  the  air  to  make  good  coupling  between  the  gelatin  solution  and  the  fat  tissue.  After 
being  cooled  in  room  temperature  for  about  15  minutes,  the  gelatin  solution  was  solidified. 
During  the  experiment,  the  transducer  also  rotationally  scanned  the  sample  from  0  to  360  degrees 
with  a  step  size  of  2.25  degrees. 

The  reconstructed  image  produced  by  our  modified  backprojection  method  is  shown  in 
Fig.  4(b),  which  agrees  with  the  original  sample  very  well.  In  particular,  the  relative  locations 
and  sizes  of  those  small  thermoacoustic  sources  are  clearly  resolved  and  perfectly  match  the 
original  ones.  Fig.  4(c)  shows  a  reconstructed  profile  (solid  curve)  at  position  x  =  21  AS  mm  of 
the  image  Fig.  4(b),  which  includes  two  gelatin  sources  with  a  distance  of  about  3  mm.  Each 
gelatin  source  has  a  distinct  profile  in  the  image.  The  boundaries  between  them  are  clearly 
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imaged.  Moreover,  the  reconstructed  profile  is  in  good  agreement  with  the  original  profile 
(dashed  curve),  which  was  a  grayscale  profile  of  the  image  Fig.  4(b).  The  half-amplitude  line 
cuts  across  the  reconstructed  profile  at  points  Bi,  Ai,  A2  and  B2,  respectively.  The  distances 
|A,B,|  =  1.72  mm  and  |A2B2|  =  1.67  mm  in  the  image  are  close  to  the  original  values  of  about 


1.80  mm  and  1.60  mm,  respectively,  which  were  measured  in  the  original  objects.  Therefore,  the 
width  of  the  profile  at  the  half-amplitude  closely  measures  its  physical  size. _ 
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Fig.  4.  (a)  Photograph  of  the  cross-section  of  a 
tissue  sample;  (b)  Reconstructed  image;  (c) 
Compare  a  line  profile  (solid  curve)  of  the 
reconstructed  image  (b)  at  x  =  27.45  mm  and 
the  corresponding  grayscale  profile  (dashed 
curve)  of  the  photograph  (a). _ 


We  here  define  a  resolving  criterion  for  estimating  spatial  resolution.  The  quarter- 
amplitude  line  cuts  across  the  profiles  at  points  Ci  and  C2,  respectively,  as  shown  in  Fig.  4(c).  If 
the  right  source  moves  to  the  position  of  the  left  one,  the  reconstructed  profile  is  equal  to  the 
spatial  summation  of  the  profiles  of  the  two  sources,  because  of  the  linear  superposition  property 
of  acoustic  waves.  When  point  C2  encounters  C\,  the  new  amplitude  at  C2  or  Ci  reaches  half 
amplitude,  and  the  two  sources  can  still  be  differentiated.  If  the  right  one  moves  more  to  the  left, 
the  new  amplitude  between  their  overlap  regions  goes  up  more  than  half  amplitude.  When  we 
use  a  half-amplitude  line  to  cut  across  the  profiles,  we  get  only  two  points  on  the  far  side  of  each 
profile,  which  means  that  these  two  sources  can  no  longer  be  clearly  distinguished.  Further, 
when  point  Ai  touches  A2,  these  two  sources  join  as  an  object. 

Therefore,  the  minimum  distance  that  can  be  differentiated  is  approximately  equal  to  the 
summation  of  the  horizontal  distance  between  point  Ai  and  Ci  and  the  horizontal  distance 
between  point  A2  and  C2.  We  have  checked  additional  pairs  of  sources  resembling  those  in  the 
image  of  Fig.  4(b),  and  found  that  this  minimum  distance  is  less  than  0.5  mm.  We  can,  therefore, 
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claim  the  spatial  resolution  in  our  experiments  reaches  less  than  0.5  mm,  which  agrees  with  the 
theoretical  spatial  resolution  limit  for  1.8  MHz  signals  whose  half  wavelength  is  ~0.5  mm  with 
the  sound  speed  of  1 .5  mm/ps. 

Of  course,  the  detecting  transducer  has  a  finite  physical  size.  If  it  is  close  to  the 
thermoacoustic  sources,  it  cannot  be  approximated  as  a  point  detector.  Its  size  will  blur  the 
images  and  decrease  the  spatial  resolution.  Therefore,  in  experiments,  the  transducer  must  be 
placed  some  distance  away  from  the  tissue  samples.  In  general,  due  to  the  finite  size  of  the 
transducer,  the  farther  away  the  transducer  is  from  the  detection  center,  the  better  the  resolution 
at  the  expense  of  the  signal. 

Other  factors  limiting  spatial  resolution  are  the  duration  of  the  microwave  pulse  and  the 
impulse  response  of  the  transducer.  In  general,  using  a  shorter  microwave  pulse  will  produce 
more  high-frequency  components  in  the  thermoacoustic  signals.  Selection  of  the  duration  of  the 
pulse  is  dependent  on  the  experimental  conditions  and  measurement  systems. 


Images  oj_ thick  samples 

The  diagram  and  the  photograph  of  a  thick  sample  are  shown  in  Fig.  5(a)  and  (b),  respectively. 
The  reconstructed  image  produced  by  our  modified  backprojection  method  is  shown  in  Fig.  5(c), 
which  agrees  with  the  original  sample  very  well.  The  relative  locations  and  sizes  of  those 
thermoacoustic  sources  perfectly  match  the  buried  objects  in  the  original  sample. 


Statement  of  Work 

Task  1:  Setting  up  the  scanning  microwave-induced-acoustic  tomography  (SMIAT)  instrument, 
Months  1-12: 


a.  Modify/connect  the  microwave  generator  and  the  ultrasonic  scanner. 

b.  Image  biological  tissues  in  vitro  with  SMIAT. 

Task  2.  Extensive  evaluation  and  optimization  of  the  SMIAT  setup,  Months  13-36: 

a.  Simulate  microwave-induced-acoustic  signals  to  provide  guidance  on  the 
experiments. 

b.  Optimize  the  ultrasonic  and  microwave  parameters  for  good  resolution  and  signal-to- 
noise  ratio. 

c.  Quantify  the  maximum  imaging  depth  with  SMIAT. 

d.  Image  biological  tissues  in  vitro  with  SMIAT  and  quantify  the  imaging  resolution. 

e.  Image  biological  tissues  in  vitro  with  SMIAT  and  ultrasonography  and  quantify  the 
contrast  improvement  of  SMIAT  over  ultrasonography. 

f.  Co-register  the  SMIAT  images  with  the  conventional  ultrasonograms. 

Task  1  has  been  successfully  accomplished  and  Task  2  is  ongoing  as  planeed.  We  are 
well  prepared  to  continue  the  project. 


Key  Research  Accomplishments 

•  An  exact  and  an  approximate  time-domain  reconstruction  algorithm  for  thermoacoustic 
tomography  in  a  spherical  geometry  was  derived  and  published. 

•  An  exact  frequency-domain  reconstruction  algorithm  for  thermoacoustic  tomography  in  a 
planar  geometry  was  derived  and  published. 

•  An  exact  frequency-domain  reconstruction  algorithm  for  thermoacoustic  tomography  in  a 
cylindrical  geometry  was  derived  and  published. 

•  High-resolution  and  high-contrast  images  were  obtained  and  published. 


Reportable  Outcomes 

Peer-reviewed  journal  articles 

1.  M.  Xu  and  L.-H.  Wang,  "Time-domain  reconstruction  for  thermoacoustic 
tomography  in  a  spherical  geometry,"  IEEE  Transactions  on  Medical  Imaging  21 
(7),  814-822  (July  2002). 

2.  Y.  Xu,  D.  Feng,  and  L.-H.  Wang,  "Exact  frequency-domain  reconstruction  for 
thermoacoustic  tomography — I:  Planar  geometry,"  IEEE  Transactions  on  Medical 
Imaging  2 1  (7),  823-828  (July  2002). 

3.  Y.  Xu,  M.  Xu,  and  L.-H.  Wang,  "Exact  frequency-domain  reconstruction  for 
thermoacoustic  tomography — II:  Cylindrical  geometry,"  IEEE  Transactions  on 
Medical  Imaging  21  (7),  829-833  (July  2002). 


13 


Conference  presentations/proceedings 

1.  M.  Xu,  X.  Wang  and  L.-H.  Wang,  "RF-  and  laser-induced  thermoacoustic 
tomography,"  pp.  406-408,  Advances  in  Optical  Imaging  and  Photon  Migration, 
Optical  Society  of  America,  Miami  Beach,  Florida,  April  7-10,  2002. 


Invited  talks 

1.  02/27/2002,  Ultrasound-mediated  biophotonic  imaging.  Laser  Diagnostic 
Technologies,  Inc.,  San  Diego,  California. 

2.  03/13/2002,  Non-invasive  skin  detection  by  oblique-incidence  reflectometry. 
STARTech  Early  Ventures,  Richardson,  Texas. 

3.  04/26/2002,  guest  speaker  for  War  On  Cancer,  7th  Annual  Relay  for  Life, 
American  Cancer  Society,  College  Station,  Texas.  The  other  guest  speakers  are 
Texas  A&M  University  System  Chancellor  Howard  Graves  and  his  wife,  Mrs. 
Gracie  Graves. 

4.  05/13/2002,  Biophotonic  imaging:  contrast,  resolution,  and  ultrasonic  mediation. 
Photonify  Technologies  Inc.,  Fremont,  California. 

5.  06/02/2002,  Ultrasound-mediated  biophotonic  imaging.  Workshop  on  Microwave 
Photonics  for  Medical  Imaging,  International  Microwave  Symposium,  Seattle, 
Washington. 

6.  06/02/2002,  Biophotonic  imaging:  contrast,  resolution,  and  ultrasonic  mediation. 
Center  for  Industrial  and  Medical  Ultrasound,  University  of  Washington,  Seattle, 
Washington. 

7.  06/1 1/2002,  Biophotonic  imaging:  contrast,  resolution,  and  ultrasonic  mediation. 
Department  of  Chemistry,  University  of  California,  Davis,  California. 

8.  06/23/2002,  Ultrasound-mediated  biophotonic  imaging.  International  Quantum 
Electronics  Conference  (IQEC)  and  the  Conference  on  Lasers,  Applications,  and 
Technologies  (LAT),  Presidium  Building  of  the  Russian  Academy  of  Sciences 
(RAS),  Moscow,  Russia. 

9.  06/27/2002,  Biophotonic  imaging:  contrast,  resolution,  and  ultrasonic  mediation. 
ESPCI,  Paris,  France. 

Degrees 

1 .  None. 
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Conclusions 

Since  October  2001,  we  have  published  three  peer-reviewed  journal  articles  in  IEEE 
Transactions  on  Medical  Imaging — a  top  imaging  journal,  published  one  conference  proceedings 
article,  and  delivered  9  invited  talks.  This  research  has  been  exciting.  We  are  well  prepared  to 
continue  the  project  with  greater  successes. 

The  combination  of  ultrasound  and  microwave  has  provided  us  a  unique  opportunity  for 
early-cancer  imaging  with  high  resolution  and  high  contrast.  We  have  made  significant  technical 
progress  in  thermoacoustic  imaging  including  data  acquisition  and  imaging  reconstruction. 
Specifically,  our  accomplishments  include  (1)  an  exact  and  an  approximate  time-domain 
reconstruction  algorithm  for  thermoacoustic  tomography  in  a  spherical  geometry  was  derived  and 
published,  (2)  an  exact  frequency-domain  reconstruction  algorithm  for  thermoacoustic 
tomography  in  a  planar  geometry  was  derived  and  published,  (3)  an  exact  frequency-domain 
reconstruction  algorithm  for  thermoacoustic  tomography  in  a  cylindrical  geometry  was  derived 
and  published,  and  (4)  high-resolution  and  high-contrast  images  were  obtained  and  published. 
The  reconstruction  is  an  inverse  source  problem  similar  to  that  in  PET  (positron  emission 
tomography);  however,  the  reconstruction  in  PET  is  based  on  geometric  optics  whereas  the 
reconstruction  in  thermoacoustic  imaging  is  based  on  diffractive/wave  optics.  We  have 
successfully  imaged  biological  tissue  with  high  resolution  and  high  contrast.  We  will  advance 
this  technology  toward  clinical  applications. 
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Time-Domain  Reconstruction  for  Thermoacoustic 
Tomography  in  a  Spherical  Geometry 

Minghua  Xu  and  Lihong  V.  Wang* 


Abstract — Reconstruction-based  microwave-induced  thermoa¬ 
coustic  tomography  in  a  spherical  configuration  is  presented. 
Thermoacoustic  waves  from  biological  tissue  samples  excited 
by  microwave  pulses  are  measured  by  a  wide-band  unfocused 
ultrasonic  transducer,  which  is  set  on  a  spherical  surface  enclosing 
the  sample.  Sufficient  data  are  acquired  from  different  directions 
to  reconstruct  the  microwave  absorption  distribution.  An  exact 
reconstruction  solution  is  derived  and  approximated  to  a  modified 
backprojection  algorithm.  Experiments  demonstrate  that  the 
reconstructed  images  agree  well  with  the  original  samples.  The 
spatial  resolution  of  the  system  reaches  0.5  mm. 

Index  Terms — Microwave,  reconstruction,  thermoacoustic, 
tomography. 

I.  Introduction 

PULSED-MICROWAVE-INDUCED  thermoacoustic  to¬ 
mography  in  biological  tissues  combines  the  advantages 
of  pure  microwave  imaging  [l]-[3]  and  pure  ultrasound 
imaging  [4],  [5].  The  wide  range  of  microwave  absorption 
coefficients  found  in  different  kinds  of  tissue  leads  to  a  high 
imaging  contrast  for  biological  tissues.  However,  it  is  difficult 
to  achieve  good  spatial  resolution  in  biological  tissues  using 
pure  microwave  imaging  because  of  the  long  wavelength  of 
microwaves.  This  problem  can  be  overcome  by  the  use  of  mi¬ 
crowave-induced  thermoacoustic  waves.  Because  the  velocity 
of  acoustic  waves  in  soft  tissue  is  ~  1 .5  mm//is,  thermoacoustic 
signals  at  megahertz  can  provide  millimeter  or  better  spatial 
resolution. 

In  thermoacoustic  tomography,  a  short-pulsed  microwave 
source  is  used  to  irradiate  the  tissue.  The  relatively  long 
wavelength  of  the  microwave,  e.g.,  ~3  cm  at  3  GHz  in  tissues, 
serves  to  illuminate  the  tissue  homogeneously.  A  wide-band 
ultrasonic  transducer  can  then  be  employed  to  acquire  the 
thermoacoustic  signals  excited  by  thermoelastic  expansion, 
which  carries  the  microwave  absorption  property  of  the  tissue. 
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The  ultrasonic  transducer  is  very  sensitive  in  detecting  small 
thermoacoustic  vibrations  from  an  object. 

The  key  problem  with  this  technique  is  how  to  determine 
the  microwave  absorption  distribution  from  the  measured 
data,  i.e.,  how  to  map  the  inhomogeneity  of  the  tissue.  One 
approach  is  to  use  focused  ultrasonic  transducers  to  localize 
the  thermoacoustic  sources  in  linear  or  sector  scans  and  then 
construct  the  images  directly  from  the  data  as  is  often  done  in 
pulse-echo  ultrasonography  [6],  [7].  An  alternative  method  is 
to  use  wide-band  point  detectors  to  acquire  thermoacoustic  data 
and  then  reconstruct  the  microwave  absorption  distribution. 
To  date,  we  have  not  seen  an  exact  inverse  solution  for  this 
specific  problem,  although  some  researchers  have  arrived  at 
approximate  reconstruction  algorithms,  such  as  the  weighted 
delay-and-sum  method  [8],  the  optimal  statistical  approach  [9], 
and  other  approach  [10]. 

Based  on  spherical  harmonic  functions,  in  this  paper  we  first 
deduce  an  exact  solution  to  the  problem  in  three-dimensional 
spherical  geometry,  which  can  be  carried  out  in  the  frequency 
domain  [  1 1  ]-[  1 4] .  The  exact  reconstruction  algorithms  in  planar 
and  cylindrical  geometries  are  reported  in  the  companion  pa¬ 
pers  [15],  [16].  Spherical  measurement  geometry  may  be  more 
suitable  for  investigation  of  external  organs  such  as  the  breast. 
We  assume  that  the  wide-band  unfocused  ultrasonic  transducer 
is  set  on  a  spherical  surface,  which  encloses  the  sample  under 
investigation.  The  data  acquired  from  different  directions  are 
sufficient  to  allow  us  to  reconstruct  the  microwave  absorption 
distribution. 

In  many  cases,  the  diameter  of  the  sphere  of  detection  is  much 
larger  than  the  ultrasonic  wavelength.  As  a  result,  an  approximate 
algorithm  can  be  deduced,  which  is  a  modified  backprojection  of 
a  quantity  related  to  the  thermoacoustic  pressure.  This  approxi¬ 
mate  algorithm  can  be  carried  out  in  the  time  domain  and  is  much 
faster  than  the  exact  solution.  In  our  initial  investigations,  we  have 
also  tested  tissue  samples  in  a  circular  measurement  configura¬ 
tion.  These  experiments  demonstrate  that  the  images  calculated 
by  the  modified  backprojection  method  agree  well  with  the  orig¬ 
inal  samples.  Moreover,  the  images  have  both  the  high  contrast 
associated  with  pure  microwave  imaging  and  the  0.5-mm  spatial 
resolution  associated  with  pure  ultrasonic  imaging. 

II.  THEORY 

A.  Fundamental  of  Thermoacoustics 

Thermoacoustic  theory  has  been  discussed  in  many  literature 
reviews  such  as  [13].  Here,  we  briefly  review  only  the  funda¬ 
mental  equations.  If  the  microwave  pumping  pulse  duration  is 
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much  shorter  than  the  thermal  diffusion  time,  thermal  diffusion 
can  be  neglected;  consequently,  the  thermal  equation  becomes 

PCp^T(r,t)  =  H(r,t)  (1) 

where  p  is  the  density;  Cp  is  the  specific  heat;  T( r,  t)  is  the 
temperature  rise  due  to  the  energy  pumping  pulse;  and  H( r,  t) 
is  the  heating  function  defined  as  the  thermal  energy  per  time 
and  volume  deposited  by  the  energy  source.  We  are  initially 
interested  in  tissue  with  inhomogeneous  microwave  absorption 
but  a  relatively  homogeneous  acoustic  property.  The  two  basic 
acoustic  generation  equations  in  an  acoustically  homogeneous 
medium  are  the  linear  inviscid  force  equation 


d2 

P  Qp  u(r,  t)  =  —  Vp(r,  t ) 


and  the  expansion  equation 


V  •  u(r,  t)  =  -^p-+pT{ r,  t) 
pc 


(2) 


(3) 


where  /?  is  the  isobaric  volume  expansion  coefficient;  c  is  the 
sound  speed;  u(r,  t)  is  the  acoustic  displacement;  and  p(r,  t) 
is  the  acoustic  pressure. 

Combining  (1)— (3),  the  pressure  p(r,  t)  produced  by  the  heat 
source  H{ r,  t)  obeys  the  following  equation: 


V2p(r,  t) 


1  d2 


<4) 


The  solution  based  on  Green’s  function  can  be  found  in  the  lit¬ 
erature  of  physics  or  mathematics  [12],  [14].  A  general  form  can 
be  expressed  as 

/? 


P( r,  = 


4tt  C„ 


dV 

dH( r',  t') 

Ir-r'l 

dtf 

|r-r'|/c) 


(5) 

The  heating  function  can  be  written  as  the  product  of  a  spatial 
absorption  function  and  a  temporal  illumination  function 

H{ r,  t)  =  (6) 

Thus,  p( r,  f)can  be  expressed  as 

^•t)=^oJS  V^r\A{r')nt)  m 

where  /'(£')  =  dl{tf)/dtf. 


B.  Exact  Reconstruction  Theory 

We  first  solve  the  problem  where  the  pulse  pumping  is  a  Dirac 
delta  function 


/(*)=*(*). 


(8) 


Suppose  the  detection  point  on  the  spherical  surface  r  =  r0, 
which  encloses  the  sample  (Fig.  1).  By  dropping  the  primes,  (7) 
may  be  rewritten  as 


p(r0,  t)  =  r)  IIJd3rAM 


A  -  ,  (9) 

47r|r0  -  r| 

where  r)  =  /i/Cp.  The  inverse  problem  is  to  reconstruct  the  ab¬ 
sorption  distribution  A(r)  from  a  set  of  datap(ro,  t)  measured 
at  positions  r0.  Taking  the  Fourier  transform  on  variable  t  of  (9), 
and  denoting  k  —  w/c,  we  get 

p(r0,  u>)  =  -iun,  JJJ  dsrA(r) 


(10) 


Fig.  1.  Acoustic  detection  scheme.  The  ultrasonic  transducer  at  position  r0 
records  the  thermoacoustic  signals  on  a  spherical  surface  with  radius  |r  —  r0|. 


where  the  following  Fourier  transform  pair  exists: 

/+oo 

p(r0,  t)exp(iut)dt, 


(11a) 


2  r-r  oo 

P(ro,t)=:7r  /  p(r0,w)exp(-Mt)du>.  (lib) 
2tt  y_o o 

We  next  derive  the  exact  solution  using  the  spherical  har¬ 
monic  function  basis.  In  the  derivation,  we  referred  to  the 
mathematical  techniques  for  ultrasonic  reflectivity  imaging 
[11].  The  mathematics  utilized  can  also  be  found  routinely 
in  the  mathematical  literature,  such  as  [12].  Here,  we  list  the 
identities  (12a}-(12f)  used  in  the  subsequent  deduction: 

1 )  The  complete  orthogonal  integral  of  spherical  harmonics 
*T(0O,  tpo) 


II 


YT(0 o,  ¥>o)nn>o,  <Po)dOo  =  Si,k6m,9 


(12a) 


where  dtt o  =  sin  Oq  dOo  d<p o  and  *  denotes  the  complex 
conjugate. 

2)  The  Legendre  polynomial 

A>j r  *^*f.  ^ 

Pi(n  •  n0)  =  Y,  Yr\e,  <p)Ytm  (0O,  <Po)  (12b) 

m=— / 

where  the  unit  vectors  n  and  no  point  in  the  directions 
( 0 ,  ip)  and  (0O,  <po)*  respectively. 

3)  The  orthogonal  integral  of  Legendre  polynomials,  derived 
from  (12a)  and  (12b) 

Jfdn0Pt( n  ■  n0)Pm(n'  •  n0)  =  SlrnP,(n  ■  n')  (12c) 

where  the  unit  vector  n'  =  r'/r'  points  in  the  direction 

(V,  h?Y 

4)  The  expansion  identity 

- 1  £( 

(*  >  0)  (12d) 

where  n  =  r/r,  no  =  r 0/r0,  ji(-)  and  /ij1^*)  are  the 


spherical  Bessel  and  Hankel  functions,  respectively. 
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5)  The  complete  orthogonal  integral  of  Bessel  functions 

A  +  OO 

J  dkk2jm(kr)jm(kr')  =  -S(r-  r’).  (12e) 

6)  The  summation  identity  of  Legendre  polynomials 

£  (2m  +  l)Pm(n  ■  n')  =  ^  ^  ^ .  (12f) 


m=0 


sin# 


First,  substituting  (1 2d)  into  (10),  we  obtain 

p(r0,  «)  =  ^  JJJ d\A( r)  £  (2 1  +  1  )j,(kr) 

■h\1]  (kro)Pi(n  •  n0).  (13) 

Then,  multiplying  both  sides  of  (13)  by  Pm(n/  •  no),  and  inte¬ 
grating  with  respect  to  no  over  the  surface  of  the  sphere,  and 
considering  the  identity  (12c),  we  obtain 

JJ  d%p( r0,  oj)Pm(ri  •  n0) 

fff  d?rA( r)  (21  +  1  )ji(kr)h\1)  (kr0) 

JJJ  l=0 

JJ dQ0Pi(n  ■  n0) Pm(ri  •  n0) 

!  JJJ d\A(r)  £(2l+l)jl(kr)h\i\krQ)^i 

■  6tmPi(n  •  n') 

=  k2t]c  JJJ d3rA(r)jm(kr)h$(kr0)Pm(n  ■  n') 


U>kl] 

47T 


wkr] 

4tc 


i.e., 


4nS(0  -  6')6(<p  -  <gQ  7T  _  , 
o„.2  ™  7  7  > 


sin0 


2  r- 


This  is  the  exact  inverse  solution  of  (9).  It  involves  summation 
of  a  series  and  may  take  much  time  to  compute.  Therefore,  it  is 
desirable  to  further  simplify  the  solution. 

C.  Modified  Backprojection 

In  experiments,  the  detection  radius  r’0  is  usually  much  larger 
than  the  wavelengths  of  the  thermoacoustic  waves  that  are  useful 
for  imaging.  Because  the  low-frequency  component  of  the  ther¬ 
moacoustic  signal  does  not  significantly  contribute  to  the  spatial 
resolution,  it  can  be  removed  by  a  filter.  Therefore,  we  can  as¬ 
sume  |A;|ro  >  1  and  use  the  asymptotic  form  of  the  Hankel 
function  to  simplify  (15).  The  following  two  identities  are  in¬ 
volved  [12]: 

1)  The  expansion  identity  similar  to  (12d) 

4?r|ro  -  r|  4tt  ^ 

■h%>(kr0)Pm( n-no),  (*  >  0).  (16a) 

2)  The  approximation  when  |fc|r0  >  1 

A£)(*ro)  *  iwof + °  (w))  (16b) 

where  h\2\-)  is  the  spherical  Hankel  function  of  the 
second  kind. 

Substituting  (16b)  into  (15),  we  get 

A(r)  «  2^^.  [[  d^°  JQ  dkP(r 0’  w)*2?o  ( 2m  + 


Qq 


[[  dfl0p( ro,  oj)Pm( n'  •  n0)  (1/ 

JJ  hrn  (kv o) 

=  k2r)C  JJJ  d3rA(r)jm  (kr)Pm  (n  •  n' ) .  (14) 

Further,  multiplying  both  sides  of  (14)  by  jm(kr '),  integrating 
them  with  respect  to  k  from  zero  to  -boo,  and  then  multiplying 
both  sides  of  (14)  again  by  (2m  + 1)  and  summing  m  from  zero 
to  oo,  and  considering  the  identity  (12e)  and  (12f),  we  get 

[[ dSl„  /+“  «,)  £  (2m+^"(*r')  P„(»'  ■  no) 

JJ  Jo  m=0  (^ro) 


■jm(kr)h^>(kr0)Pm(n  ■  n0).  (17) 

Considering  the  form  of  (16a),  the  above  equation  can  be 
rewritten  as 


Mr)  = 


»o 


27T27]C 


#p+oo 

dtooj  < 

Q0 


dkp( r0,  u)(ik) 


exp(-zfc|ro  -  r|) 
|ro-r| 


S7q 


exp 


err  w  p-roo 

=  rjc  JJJ  dzrA(r)  ^  (2m  +  l)Pm(n  *  n')  J  dkk‘ 

•  jm(kr')jm(ki') 

=  rjc  JJJ  d?rA(  r) 

=  27r27/cyl(r'). 

Finally,  dropping  the  primes,  we  can  rewrite  the  equation  as 

A<-T)=*bJJdU°l  dkfin’w) 

n-no).  (15) 


ko  -  r| 

Because p(r,  t)  is  a  real  function, p*(r,  a;)  =  p(r,  -a;).  Taking 
the  summation  of  the  above  equation  with  its  complex  conjugate 
and  then  dividing  it  by  two,  we  get 

■A(r)  =  4^2vc  JJ  d^°  J  dkp(r0,  u)(ik) 

Qo 

exp(— iA:|r0  -  r|) 


|r0  • 


27T7/C3 


exp 


m=0 


h$(kr0) 


Ihrj: 

Qq 

(ggg1) 


ko  -  r| 


<Ljp( ro,  u>)(—iw) 
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Recalling  the  inverse  Fourier  transform  (lib),  we  get 


M  r)  =  - 


2tCT)(? 


// 

tto 


dil  o 


dp{ r0,  t) 


lro  “  rl 


dt 


(18) 


t=|r0~rj/c 


i.e., 


Mr)  =  - 


27T7/C4 


1  dp( r0,  t ) 


dt 


Q0 


(19) 


t=|r0— r|/c 


Equation  (19)  shows  that  the  absorption  distribution  can  be 
calculated  in  the  time  domain  by  the  means  of  backprojection 
and  coherent  summation  over  spherical  surfaces  of  the  quantity 
-(l/f)(<9p(r0,  t)/dt)  instead  of  the  acoustic  pressure  itself. 
This  approximate  algorithm  requires  less  computing  time  than 
the  exact  solution  (15). 

For  initial  investigations,  we  measure  the  samples  in  a  circular 
configuration.  In  these  cases,  the  backprojection  is  carried  out 
in  a  circle  around  the  slices,  and  (19)  can  be  simplified  to 


^l(r)  =  - 


27t?7c^ 


/ 


j  1  #p(ro,  t) 
dipo  t  dt 


(20) 


t=r|r0—  rj/c 


III.  Experimental  Method 

A.  Diagram  of  Setup 

Fig.  2  shows  the  experimental  setup  for  the  circular  measure¬ 
ment  configuration,  which  is  modified  from  our  previous  paper 
[7].  For  the  convenience  of  the  reader,  the  system  is  briefly  de¬ 
scribed  here.  The  unfocused  transducer  (V323,  Panametrics)  has 
a  central  frequency  of  2.25  MHz  and  a  diameter  of  6  mm.  It  is 
fixed  and  it  points  horizontally  to  the  center  of  the  rotation  stage, 
which  is  used  to  hold  the  samples.  For  good  coupling  of  acoustic 
waves,  both  the  transducer  and  the  sample  are  immersed  in  min¬ 
eral  oil  in  a  container. 

The  microwave  pulses  are  transmitted  from  a  3-GHz  mi¬ 
crowave  generator  with  a  pulse  energy  of  10  mJ  and  a  width 
of  0.5  /is,  and  then  delivered  to  the  sample  from  the  bottom 
by  a  rectangular  waveguide  with  a  cross  section  of  72  mm  x 
34  mm.  A  function  generator  (Protek,  B-180)  is  used  to  trigger 
the  microwave  generator,  control  its  pulse  repetition  frequency, 
and  synchronize  the  oscilloscope  sampling.  The  signal  from 
the  transducer  is  first  amplified  through  a  pulse  amplifier, 
then  recorded  and  averaged  200  times  by  an  oscilloscope 
(TDS640A,  Tektronix).  A  personal  conputer  is  used  to  control 
the  step  motor  for  rotating  the  sample  and  transferring  the  data. 

Last,  we  want  to  point  out  that,  in  our  experiments,  the 
smallest  distance  r0  between  the  rotation  center  and  the 
surface  of  the  transducer  is  4.3  cm.  In  the  frequency  domain 
(100  KHz- 1.8  MHz),  |fc|r0  =  27rr0//cwith  1 .5  mm/^s,  we  get 
18  <  |fc|r0  <  330.  Therefore,  the  required  condition  \k\r0  >  1 
for  the  modified  backprojection  algorithm  is  satisfied. 

B.  Technical  Consideration 

During  measurement,  we  find  that  the  piezoelectric  signal 
50(r0,  t)  detected  by  the  transducer  includes  the  thermal 
acoustic  signal  ^(ro,  t)  as  well  as  some  noise.  The  noise 
comes  from  two  contributors.  One  is  the  background  random 
noise  of  the  measurement  system,  which  can  be  suppressed  by 
averaging  the  measured  data.  The  other  part,  5mp(f),  results 
from  the  microwave  pumping  via  electromagnetic  induction. 


Fig.  3.  (a)  The  temporal  profile  of  the  microwave  pulse;  (b)  the  temporal 
profile  of  the  impulse  response  of  the  transducer;  (c)  compare  the  normalized 
amplitudes  of  the  spectrum  /(/)/?(/),  G(f)  and  fG(f). 
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Fig.  4.  (a)  An  example  of  temporal  piezoelectric  signal;  (b)  an  example  of 
temporal  noise;  (c)  an  example  of  filtered  signal;  and  (d)  an  example  of  filtered 
thermoacoustic  signals  detected  at  different  angular  positions  from  0°  to  360°. 


The  pumping  component  of  the  noise  can  be  measured  without 
a  sample,  and  can  be  subtracted  from  the  measured  data 

S(r0,  t)  w  50(r0,  t)  -  Smp{t).  (21) 

In  fact,  the  transducer  is  not  a  real  point  detector.  For  sim¬ 
plicity,  we  can  ignore  its  size  if  we  put  it  far  away  from  the 
sample.  However,  we  still  have  to  consider  the  impulse  response 
R(t)  of  the  transducer  and  the  pumping  duration  I(t )  of  the  mi¬ 
crowave  pulse.  In  general,  the  measured  thermoacoustic  signal 
can  be  written  as  a  convolution 

S(r0,  t)  =  p(r0,  t)  *  1(f)  *  R(t)  (22) 

where  p(ro,  f)  is  the  thermoacoustic  signal  with  delta-pulse  mi¬ 
crowave  pumping.  In  the  frequency  domain,  (22)  can  be  written 
as 

5(r0,  u>)  =  p(r0,  u>)I(w)R(u>)  (23) 

where 


r f«> 

L 

I(t)exp(ivt)  dt , 

(24a) 

R{")  = 

L 

R(t)  exp(zo;f)  dt . 

(24b) 

Because  of  the  presence  of  noise  and  the  finite  bandwidth  of 
I(w)  and  R(u),  an  appropriate  deconvolution  algorithm  should 
be  used  to  calculate  p(r0,  tu).  In  the  reconstruction,  only  the 
high-frequency  component  of  the  thermoacoustic  signal  is  re¬ 
quired.  Therefore,  we  compute  p(ro,  u>)Cx(u>)  instead,  where 
G(u)  is  a  high-frequency  bandpass  filter  such  as  a  Gaussian 
filter 

G(aO=exp[— ] 

and  a  and  w o  are  two  parameters  of  the  filter,  =  27r/  and 
u0  =  2? r/0. 

In  our  experiments,  1(f)  is  approximately  a  rectangular  func¬ 
tion  with  duration  r  =  0.5  /is  and  its  temporal  profile  is  shown 
in  Fig.  3(a).  Its  spectrum  I(u)  covers  the  range  from  0  to  2  MHz. 
The  transducer  that  we  used  is  of  the  videoscan  type  with  a  cen¬ 
tral  frequency  of  fo  =  2.25  MHz,  and  the  temporal  profile  of 
the  impulse  response  is  shown  in  Fig.  3(b).  It  is  observed  that 
the  generated  thermoacoustic  signal  under  microwave  pumping 
with  duration  r  =  0.5  /is  exists  primarily  in  a  frequency  range 
below  1.8  MHz.  We  chose  the  parameters  a  =  3.6  and  /0  = 
1.25  MHz  in  the  Gaussian  filter 

G(/)  =  exp[-«(i-l)] 

to  eliminate  the  noise  at  high  as  well  as  low  frequencies.  The 
spectrum  G(f)  is  shown  as  the  dash-dot  line  in  Fig.  3(c). 
We  compared  the  normalized  spectrum  I(f)R(f)  [solid  line 
in  Fig.  3(c)]  with  fG(f)  [dash  line  in  Fig.  3(c)],  and  found 
\fG(f)\  «  \I(f)R(f)\  when  /  <  2  MHz.  Of  course,  this 
approximated  equality  is  a  special  case  for  our  measurement 
system  only.  Therefore,  the  filtered  dp(r0,  t)/dt  can  be  simply 
calculated  by  an  inverse  fast  Fourier  transform  (IFFT) 

dvijZ'  » IFFT  {S(r0,  w)F(a;)}  (25) 
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Fig.  5.  (a)  Cross  section  of  a  tissue  sample;  (b)  reconstructed  image;  and 
(c)  a  line  profile  of  the  reconstructed  image  at  y  =  31.5  mm. 

where  F(w)  is  a  wide  bandpass  filter,  which  is  used  to  further 
eliminate  noise  at  high  and  low  frequencies  in  order  to  guarantee 


the  condition  |fc|ro  >■  1  for  the  modified  backprojection.  A 
simple  filter  is 


0.1  MHz  <  /  <  1.8  MHz 
otherwise. 


(26) 


IV.  Results  and  Discussion 

Finally,  we  use  the  above  modified  backprojection  algorithm 
and  the  experimental  method  to  investigate  some  tissue  samples. 

A.  Experimental  Data  Preprocessing 

The  measured  piezoelectric  data  include  the  useful  thermoa¬ 
coustic  signal  as  well  as  some  noise  data  as  illustrated  by  the  fol¬ 
lowing  example.  Fig.  4(a)  is  a  typical  measured  temporal  piezo¬ 
electric  signal,  which  is  from  the  sample  shown  in  Fig.  5(a). 
One  portion  of  the  noise  resulting  from  the  microwave  pumping 
looks  like  the  curve  in  Fig.  4(b),  which  is  acquired  at  the  same 
sampling  rate  and  the  same  delay  time  with  the  transducer  in 
the  same  position  as  the  curve  in  Fig.  4(a).  Because  the  slice  is 
very  thin,  the  thermoacoustic  signal  is  not  much  higher  than  the 
noise  resulting  from  the  microwave  pumping.  Next,  we  subtract 
the  noise  from  the  raw  thermoacoustic  signal  and  use  a  wide 
bandpass  filter  to  eliminate  some  of  the  useless  low-frequency 
and  high-frequency  components.  This  processed  data  is  shown 
in  Fig.  4(c);  it  is  much  cleaner  than  the  raw  data  in  Fig.  4(a).  The 
filtered  thermoacoustic  signals  detected  at  different  angular  po¬ 
sitions  from  0°  to  360°  are  shown  in  Fig.  4(d). 

B.  Image  Contrast 

Image  contrast  is  an  important  index  for  biological  imaging. 
Fig.  5(a)  shows  a  tested  sample,  which  was  photographed  after 
the  experiment.  The  sample  was  made  according  to  the  fol¬ 
lowing  procedure.  First,  we  cut  a  thin  piece  of  homogeneous 
pork  fat  tissue  and  shaped  it  arbitrarily  to  form  a  base.  Its  thick¬ 
ness  is  5  mm  and  its  maximum  diameter  is  4  cm.  Then  we  used 
different  screwdrivers  to  carefully  make  two  pairs  of  holes  that 
were  approximately  4  and  6  mm  in  diameter,  respectively.  Fi¬ 
nally,  one  big  and  one  small  hole  on  the  left  side  were  filled 
with  pork  muscle,  while  the  two  holes  on  the  right  side  were 
filled  with  pork  fat  of  the  same  type  as  that  which  made  up  the 
base. 

In  the  experiment,  the  transducer  rotationally  scanned  the 
sample  from  0°  to  360°  with  a  step  size  of  2.25°.  The  detec¬ 
tion  radius  7  o  was  4.3  mm.  We  used  the  160  series  of  data  as 
shown  in  Fig.  4(d)  to  calculate  the  image  by  our  modified  back- 
projection  method. 

The  reconstructed  image  is  shown  in  Fig.  5(b).  The  outline 
and  size  of  the  fat  base  as  well  as  the  sizes  and  locations  of 
the  two  muscle  pieces  are  in  good  agreement  with  the  original 
sample  in  Fig.  5(a).  Fig.  5(c)  shows  a  line  profile  for  the  small 
piece  of  muscle  in  the  image.  It  indicates  that  the  contrast  be¬ 
tween  the  fat  and  the  muscle  is  very  high.  This  high  contrast 
is  due  to  the  low  microwave  absorption  capacity  of  fat  and  the 
high  absorption  capacity  of  muscle:  at  3  GHz,  the  penetration 
depth  for  muscle  and  fat  are  1.2  and  9  cm,  respectively.  How¬ 
ever,  the  two  pieces  of  fat  are  not  visible  in  the  image  [Fig.  5(b)], 
which  means  the  minute  mechanical  discontinuity  between  the 
boundaries  of  muscle  and  fat  does  not  contribute  much  to  the 
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thermoacoustic  signal.  On  the  contrary,  discontinuity  improves 
the  strength  of  the  echo  sounds  in  pure  ultrasound  imaging. 

C.  Spatial  Resolution 

Spatial  resolution  is  another  important  index  for  biological 
imaging.  We  used  samples  with  a  set  of  small  thermoacoustic 
sources  to  test  the  resolution.  One  tested  sample  is  shown  in 
Fig.  6(a),  which  was  also  photographed  after  the  experiment  was 
completed. 

The  sample  was  made  according  to  the  following  procedure. 
First,  we  cut  a  thin  piece  of  homogeneous  pork  fat  tissue  and 
made  it  into  an  arbitrary  shape.  Its  thickness  was  5  mm  with  a 
maximum  diameter  of  4  cm.  Then  we  used  a  small  screwdriver 
to  carefully  make  a  set  of  small  holes  about  2  mm  in  diameter.  In 
the  meantime,  we  prepared  a  hot  solution  with  5%  gelatin,  0.8% 
salt,  and  a  drop  of  dark  ink  (to  improve  the  photographic  prop¬ 
erties  of  the  sample).  Next,  we  used  an  injector  to  inject  a  drop 
of  the  gelatin  solution  into  each  small  hole  and  subsequently 
blew  out  the  air  to  make  good  coupling  between  the  gelatin  so¬ 
lution  and  the  fat  tissue.  After  being  cooled  in  room  temperature 
for  about  15  min,  the  gelatin  solution  was  solidified.  During  the 
experiment,  the  transducer  also  rotationally  scanned  the  sample 
from  0°  to  360°  with  a  step  size  of  2.25°.  The  detection  radius 
ro  was  4.3  mm. 

The  reconstructed  image  produced  by  our  modified  backpro- 
jection  method  is  shown  in  Fig.  6(b);  it  also  agrees  with  the  orig¬ 
inal  sample  well.  In  particular,  the  relative  locations  and  sizes  of 
those  small  thermoacoustic  sources  are  clearly  resolved  and  per¬ 
fectly  match  the  original  ones.  Fig.  6(c)  shows  a  reconstructed 
profile  (solid  curve)  at  position  x  =  27.45  mm  of  the  image 
Fig.  6(b),  which  includes  two  gelatin  sources  with  a  distance 
of  about  3  mm.  Each  gelatin  source  has  a  distinct  profile  in  the 
image.  The  boundaries  between  them  are  clearly  imaged.  More¬ 
over,  the  reconstructed  profile  is  in  good  agreement  with  the 
original  profile  (dashed  curve),  which  was  a  grayscale  profile 
of  the  image  Fig.  6(b).  The  half-amplitude  line  cuts  across  the 
reconstructed  profile  at  points  Bi,  Ai,  A2,  and  B2,  respectively. 
The  distances  |AiBi|  =  1.72  mm  and  |A2B2|  =  1.67  mm 
in  the  image  are  close  to  the  original  values  of  about  1 .80  and 
1.60  mm,  respectively,  which  were  measured  in  the  original  ob¬ 
jects.  Therefore,  the  width  of  the  profile  at  the  half-amplitude 
closely  measures  its  physical  size. 

We  here  define  a  resolving  criterion  to  estimate  the  spatial 
resolution.  The  quarter-amplitude  line  cuts  across  the  profiles  at 
points  Ci  and  C2,  respectively,  as  shown  in  Fig.  6(c).  If  the  right 
source  moves  to  the  position  of  the  left  one,  the  reconstructed 
profile  is  equal  to  the  spatial  summation  of  the  profiles  of  the  two 
sources,  because  of  the  linear  superposition  property  of  acoustic 
waves.  When  point  C2  encounters  Ci,  the  new  amplitude  at  C2 
or  Ci  would  reach  the  half  amplitude,  and  the  two  sources  could 
still  be  differentiated.  If  the  right  one  moves  more  to  the  left,  the 
new  amplitude  between  their  overlap  regions  would  elevate  to 
more  than  the  half  amplitude.  When  we  use  a  half-amplitude 
line  to  cut  across  the  profiles,  we  get  only  two  points  on  the  far 
side  of  each  profile,  which  means  that  these  two  sources  can  no 
longer  be  clearly  distinguished.  Further,  when  point  Ai  touches 
A2,  these  two  sources  join  as  a  single  object  in  the  image. 
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Fig.  6.  (a)  Cross  section  of  a  tissue  sample;  (b)  reconstructed  image;  and 
(c)  comparison  between  a  line  profile  (solid  curve)  of  the  reconstructed  image 
(b)  at  x  =  27.45  mm  and  the  corresponding  grayscale  profile  (dashed  curve) 
of  the  original  image  (a). 
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Fig.  7.  (a)  Diagram  of  the  sample  structure  and  the  measurement;  (b)  cross 

section  of  the  tissue  sample;  and  (c)  reconstructed  image. 


Therefore,  the  minimum  distance  that  can  be  differentiated  is 
approximately  equal  to  the  summation  of  the  horizontal  distance 
between  point  A*  and  Ci  and  the  horizontal  distance  between 


point  A 2  and  C2.  We  have  checked  additional  pairs  of  sources 
resembling  those  in  the  image  of  Fig.  6(b),  and  found  that  this 
minimum  distance  is  less  than  0.5  mm.  We  can,  therefore,  claim 
that  the  spatial  resolution  in  our  system  reaches  less  than  0.5 
mm,  which  agrees  with  the  theoretical  spatial  resolution  limit 
for  1 .8-MHz  signals  whose  half  wavelength  is  ~0.5  mm  with 
the  sound  speed  of  1 .5  mm///s. 

Of  course,  the  detecting  transducer  has  a  finite  physical  size. 
If  it  is  close  to  the  thermoacoustic  sources,  it  cannot  be  approx¬ 
imated  as  a  point  detector.  Its  size  will  blur  the  images  and 
decrease  the  spatial  resolution.  Therefore,  in  experiments,  the 
transducer  must  be  placed  some  distance  away  from  the  tissue 
samples.  In  general,  due  to  the  finite  size  of  the  transducer,  the 
farther  away  the  transducer  is  from  the  detection  center,  the 
better  the  resolution  at  the  expense  of  the  signal  strength. 

Other  limiting  factors  of  spatial  resolution  include  the 
duration  of  the  microwave  pulse  and  the  impulse  response  of 
the  transducer.  In  general,  using  a  shorter  microwave  pulse 
will  produce  more  high-frequency  components  in  the  thermoa¬ 
coustic  signals.  The  disadvantages  resulting  from  employing  a 
shorter  pulse,  however,  are  insufficient  energy  delivery  and  a 
decrease  in  the  signal-to-noise  ratio.  Selection  of  the  duration 
of  the  pulse  is  dependent  on  the  experimental  conditions 
and  measurement  systems.  In  biological  tissues,  microwaves 
at  300  MHz~3  GHz  with  0. 1~1  ^s  pulse  width  are  often 
adopted.  Therefore,  the  high-frequency  of  the  thermoacoustic 
signals  reaches  several  MHz.  Such  a  wide-band  transducer  for 
measuring  acoustic  waves  at  ~MHz  is  widely  available. 

D.  Thick  Sample 

The  advantage  of  using  microwave  is  its  long  penetration 
depth  in  soft  tissue.  A  microwave  can  reach  a  tumor  buried  in¬ 
side  tissue  and  heat  it  to  generate  thermoacoustic  waves.  One 
tested  sample  is  shown  in  Fig.  7(a).  The  experiment  was  con¬ 
ducted  according  to  a  procedure  similar  to  the  one  above.  Three 
small  absorbers  were  buried  inside  a  big  fat  base.  The  big  pork 
fat  tissue  had  a  maximum  diameter  of  7  cm.  Screwdrivers  were 
used  to  carefully  make  three  holes  about  5  mm  in  diameter  with 
a  depth  of  2.5  cm.  Next,  an  injector  was  used  to  inject  a  drop  of 
the  same  gelatin  solution  as  above  into  each  small  hole,  and,  sub¬ 
sequently,  air  was  blown  out  to  improve  the  coupling  between 
the  gelatin  solution  and  the  fat  tissue.  These  gelatin  sources  were 
about  5  mm  in  diameter.  After  being  cooled  at  room  tempera¬ 
ture  for  about  15  min,  the  gelatin  solutions  solidified.  The  pho¬ 
tograph  of  the  sample  at  this  stage  is  shown  in  Fig.  7(b).  Finally, 
the  holes  were  filled  with  fat,  and  the  gelatin  sources  were  buried 
in  the  fat  tissue. 

During  the  experiment,  a  microwave  was  transmitted  out  to 
the  sample  from  below.  The  transducer  rotationally  scanned  the 
sample,  including  the  gelatin  sources,  from  0°  to  360°  in  a  plane 
as  Fig.  7(a)  shows.  The  distance  between  the  transducer  and  the 
rotation  center  was  7  cm.  The  reconstructed  image  produced  by 
our  modified  backprojection  method,  which  agrees  well  with 
the  original  sample,  is  shown  in  Fig.  7(c). 

The  above  experiments  verified  the  principle  of  the  modi¬ 
fied  backprojection  algorithm,  which  implies  back  projection 
and  coherent  summation  over  spherical  surfaces.  In  particular,  a 
set  of  circular  measurement  data  would  be  sufficient  to  yield  a 
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satisfactory  cross-sectional  image  for  a  sample  with  only  small 
absorption  sources  in  the  same  horizontal  plane  and  a  lower 
absorption  background.  Of  course,  for  a  complicated  sample, 
data  from  only  a  circular  measurement  would  be  insufficient 
for  3-D  reconstruction  unless  cylindrical  focusing  is  employed. 
This  limited  view  problem  will  be  addressed  in  our  future  work. 

Finally,  we  must  point  out  that  an  inhomogeneous  acoustic 
property,  such  as  the  speed  of  sound  variation,  might  result  in  re¬ 
construction  errors.  Fortunately,  the  speed  of  sound  in  most  soft 
tissue  is  relatively  constant  at  ~  1.5  mm///s.  The  above  experi¬ 
ments  demonstrated  that  the  small  speed  variations  between  fat 
and  muscle  or  gelatin  did  not  result  in  significant  reconstruction 
artifacts.  The  reason  is  that  thermoacoustic  waves  are  produced 
internally  by  microwave  absorption  and  are  propagated  one-way 
to  the  detectors.  Thus,  a  small  speed  variation  does  not  affect  the 
travel  time  of  the  sound  very  much  in  a  finite-length  path,  for  ex¬ 
ample,  10  cm,  which  is  comparable  to  a  typical  breast  diameter. 
Therefore,  in  thermoacoustic  tomography,  satisfactory  contrast 
and  resolution  are  obtainable  even  in  tissue  with  a  small  degree 
of  acoustic  inhomogeneity. 

V.  Conclusion 

Pulsed-microwave-induced  thermoacoustic  tomography  of 
inhomogeneous  tissues  has  been  studied.  Both  an  exact  inverse 
solution  and  a  modified  backprojection  algorithm  have  been 
derived,  which  are  based  on  the  data  acquired  by  wide-band 
point  detectors  on  a  spherical  surface  that  encloses  the  sample 
under  study.  A  set  of  experiments  on  tissue  samples  has  been 
investigated  under  a  circular  measurement  configuration.  The 
reconstructed  images  calculated  by  the  modified  backprojec¬ 
tion  method  agree  well  with  the  original  ones.  Results  indicate 
that  this  technique  using  reconstruction  theory  is  a  powerful 
imaging  method  that  results  in  good  contrast  and  good  spatial 
resolution  (0.5  mm),  which  can  be  used  for  the  investigation  of 
tissues  with  inhomogeneous  microwave  absorptions. 
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Exact  Frequency-Domain  Reconstruction  for 
Thermoacoustic  Tomography — I:  Planar  Geometry 
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Abstract — We  report  an  exact  and  fast  Fourier-domain  re¬ 
construction  algorithm  for  thermoacoustic  tomography  in  a 
planar  configuration  assuming  thermal  confinement  and  constant 
acoustic  speed.  The  effects  of  the  finite  size  of  the  detector  and 
the  finite  length  of  the  excitation  pulse  are  explicitly  included  in 
the  reconstruction  algorithm.  The  algorithm  is  numerically  and 
experimentally  verified.  We  also  demonstrate  that  the  blurring 
caused  by  the  finite  size  of  the  detector  surface  is  the  primary 
limiting  factor  on  the  resolution  and  that  it  can  be  compensated 
for  by  deconvolution. 

Index  Terms — Fourier-domain  reconstruction,  planar,  thermo¬ 
acoustic  tomography. 

I.  Introduction 

USING  thermoacoustic  tomography  (TAT)  to  image  bio¬ 
logical  tissues  has  two  primary  advantages.  The  first  is 
the  high  spatial  resolution  comparable  with  pure  ultrasound 
imaging.  The  second  advantage  results  from  the  large  contrast 
in  microwave  absorption  that  exists  between  cancerous  tissue 
and  the  normal  tissue  [l]-[7].  Reviews  of  TAT  and  related 
works  [8]— [17]  can  be  found  in  [1 1]  and  [18]. 

Various  reconstruction  algorithms  for  TAT  [8],  [9],  [16],  [18], 
[19]  have  been  reported.  Under  the  approximation  that  the  dis¬ 
tance  between  the  detector  and  the  absorbing  object  is  much 
larger  than  the  dimension  of  the  absorbing  object,  a  three-di¬ 
mensional  (3-D)  Radon  transform  was  applied  to  reconstruct 
the  object  in  TAT  [8].  However,  the  fact  that  this  approxima¬ 
tion  may  not  always  hold  in  real-world  situations  limits  the  ap¬ 
plication  of  this  method.  Further,  the  spatial  resolutions  of  the 
imaging  system  using  this  reconstruction  method  are  limited  by 
blurs  [20]  caused  by  the  finite  size  of  the  transducer  surface, 
the  finite  width  of  the  stimulating  pulse,  and  the  finite  band¬ 
width  of  the  transducers  and  amplifiers.  Among  these  effects, 
the  blur  from  the  size  of  the  transducer  surface  is  expected  to  be 
the  largest  contributor  to  the  total  blur.  The  analysis  of  error  is 
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limited  in  numerical  simulations,  and,  hence,  no  analytical  form 
was  available  prior  to  this  work.  A  time-domain  beam-forming 
technique  was  applied  in  one  study  to  image  reconstruction  for 
the  photoacoustic  scanning  of  tissue  structures  [9].  A  weighted 
delay-and-sum  algorithm  was  used  to  account  for  the  near-field 
effect  and  to  reduce  noise.  This  algorithm  is  an  approximate  one, 
and  its  lateral  resolution  is  limited  by  the  size  of  the  detector  sur¬ 
face.  The  above  reconstructions  were  implemented  in  the  time 
domain  and  consequently  are  time-consuming,  especially  in  3-D 
tomography.  TAT  was  also  obtained  in  a  way  similar  to  that  used 
in  conventional  B-scan  ultrasonic  imaging,  but  it  had  difficulty 
detecting  the  boundaries  of  objects  that  are  oblique  to  the  trans¬ 
ducer  axis  [16].  Exact  reconstructions  have  been  implemented 
for  TAT  in  spherical  and  cylindrical  configurations  in  the  com¬ 
panion  papers  [18],  [19]. 

Next,  we  present  our  studies  on  an  exact  and  fast  reconstruc¬ 
tion  algorithm  using  a  Fourier  transform  for  TAT  in  a  planar 
configuration.  The  reconstruction  of  an  image  by  Fourier  trans¬ 
form  has  been  used  in  X-ray  computed  tomography  [21],  ul¬ 
trasonic  reflectivity  imaging  [22]-[24],  and  diffraction  tomog¬ 
raphy  [25]  successfully.  The  computation  complexity  is  reduced 
greatly  due  to  the  efficiency  of  the  Fourier  transform.  We  devel¬ 
oped  image  reconstruction  by  Fourier  transform  for  planar  TAT 
and  obtained  an  exact  reconstruction  algorithm  for  the  first  time. 
Furthermore,  some  limitations  from  experiments,  such  as  the  ef¬ 
fects  of  the  finite  size  of  the  detectors  and  the  finite  length  of  the 
excitation  pulse,  are  included  explicitly  in  the  reconstruction  al¬ 
gorithm.  The  reconstruction  algorithm  is  verified  by  both  nu¬ 
merically  simulated  and  experimental  results.  Our  simulations 
also  demonstrate  that  the  blur  due  to  the  finite  size  of  the  de¬ 
tector  surface,  which  is  a  key  limiting  factor  on  the  resolution 
of  images  [9],  [20],  can  be  alleviated  by  deconvolution  with  re¬ 
spect  to  the  size  of  the  detector  surface.  Other  effects  that  may 
cause  blurring  of  images  can  be  treated  in  a  similar  way.  In  our 
initial  experiments,  an  image  in  good  agreement  with  the  real 
objects  was  reconstructed  and  the  deconvolution  improved  the 
resolution  of  the  imaging  system. 

II.  Methods 
A.  Image  Reconstruction 

Assume  that  the  detector  scans  within  the  plane  2  =  0  and 
that  the  object  is  distributed  only  in  the  half  space  z'  >  0. 
In  order  to  obtain  a  spatial  resolution  of  about  1  mm,  the  mi¬ 
crowave  pulse  should  be  set  to  less  than  ~1  /is  because  the 
speed  of  sound  in  soft  biological  tissue  is  ~1.5  mm//xs.  For 
these  parameters,  the  diffusion  term  in  the  heat  conduction  equa¬ 
tion  is  about  six  orders  of  magnitude  less  than  the  term  of  the 
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first-derivative  of  the  temperature  [26].  Therefore,  heat  conduc¬ 
tion  can  be  ignored.  This  is  known  as  the  assumption  of  thermal 
confinement.  In  this  case,  the  acoustic  wave  p( r,  t)  is  related  to 
microwave  absorption  H (r,  t)  by  the  following  wave  equation 
[26]: 


2%2-vyr,  = 

C  ot 


(1) 


where  t  =  tvs,  vs  is  the  acoustic  speed,  C  is  the  specific  heat, 
and  /?  is  the  coefficient  of  volume  thermal  expansion.  In  (1),  the 
acoustic  speed  is  assumed  constant,  which  will  be  further  ad¬ 
dressed  in  the  discussion  section.  Equation  (1)  can  be  rewritten 
in  terms  of  H( r',  t ): 


P( r,  i)  = 


0v» 
4t tC 


dH(r'7 1')  dr ' 
dtf  |r  —  r'| 


(2) 


where  if  =  t  -  |r  -  r'|.  The  source  term  H(rf,  t)  can  further  be 
written  as  the  product  of  a  purely  spatial  and  a  purely  temporal 
component,  i.e., 


H( r',  t)  =  hy{r')r)(t)  (3) 


where  Jo  is  a  scaling  factor  proportional  to  the  incident  radiation 
intensity,  ^(r')  describes  the  microwave  absorption  properties 
of  the  medium  at  r'.  r/(?)  describes  the  shape  of  the  irradiating 
pulse  and  is  a  nonnegative  function  whose  integration  over  time 
equals  the  pulse  energy.  Substituting  (3)  into  (2)  results  in 


P(r,  <)  = 


lofty* 
4  nC 


dr' 

|r  —  rT 


(4) 


We  proceed  by  transforming  the  time-dependent  wave  equa¬ 
tion  into  the  temporal-frequency  domain.  Denoting  the  Fourier 
transforms  of  p  and  r/  by  p  and  r} ,  we  have 

/oo 

p(r,  k )  exp(ifct)  dk  (5) 

-OO 

/OO 

rj(k)  exp (ikt)  dk.  (6) 

-oo 


Substituting  (5)and  (6)  into  (4)  results  in 


p(r,  k)  = 


ipvsIokri(k) 
4t rC 


exp(— i/c|r  —  r'|) 
Ir-r'l 


dr'. 

(7) 


If  the  acoustic  signals  are  collected  along  a  line  or  in  a  plane, 
for  example  at  z  =  0,  following  the  line  of  Nortan  and  Linzer 
in  [22],  it  can  be  shown  that  for  the  case  |A;|  >  p  and  z'  >  0 


-p,  n  -  ^vsl()ky(k)sgii(k) 

{U ’  }  "  2Cy/k^? 

■  J  4>(ti,  v,  z')  exp  ^—iz'sgn(k)\/k2  —  p2^  dz  (8) 
where  p2  =  u2  +  v2,  sgn (k)  is  the  signum  function 
P(u,  v,  k)  =  ^  JJp(x>  V’  °>  k) 


•  exp  (-i(ux  +  vy))  dx  dy  (9) 


and 

*(U’  V’  Z')  =  {2 


•  exp  (— i(ux '  +  vy'))  dx*  dy*.  (10) 


Equation  (8)  can  further  be  simplified  to 

■n(3vsIokv(k)sga(k)$i(u,  v,  sgn(fc) \/k2-p2^J 
C\Jk2  -  pl 


P(u ,  v7  k) 
where 


(ID 


1  f°° 

$i(u,  v,  w)  =  —  /  4>(w,  vy  zf)exp(-~iwz')dz\  (12) 

2?r  J — oo 


The  lower  limit  of  the  above  integration  is  changed  from  0  to 
— oo  because  4>(it,  v7  zf)  = 0  when  z*  <  0.  Equation  (11)  gives 
an  exact  mapping  relation  between  the  spectrum  of  the  collected 
signals  and  the  spectrum  of  the  distribution  of  microwave  en¬ 
ergy  deposition  and  is  the  essence  of  our  reconstruction  method. 
However,  (11)  stands  only  if  the  acoustic  detector  is  a  point  de¬ 
tector.  In  practice,  the  detector  is  of  finite  size,  whose  surface 
shape  can  be  described  by  R(x ,  y).  The  signal  from  the  detector 
pd(xj  y,  t)  can  be  expressed  as  an  integral  of  the  acoustic  wave 
p(r ,  t)  over  the  detector  surface 

Pd{x,  y,t)  =  JJ  p(x',  y1,  t)R{x'  -  x,  y'  -  y)  dx'  dy'.  (13) 

a 

After  transforming  (13)  into  the  temporal-  and  spatial-frequency 
domain,  we  have 


Pd(u ,  v ,  k)  =  47 r2P(u,  v,  k)R(~u ,  -v)  (14) 


where  Pd(u,  v7  k)  is  the  temporal  and  spatial  Fourier  transform 
of  pd(x,  y ,  f),  and  R(u ,  v)  is  the  spatial  Fourier  transform  of 
R(x ,  y).  Substituting  (14)  into  (11)  results  in 


Pd(u,  v ,  k)  = 

AK*pvjQkTj{k)$&i{k)R{-u,  -v)$i  (u7  v7  sgn (k)  y/k2-p2^ 
C\Jk2—p2 

(15) 


Mapping  the  («,  v,  k)  space  into  the  (u,  v7  w)  space  by  the 
relation 

w  =  sgn  {k)\Jk2  —  p2  (16) 

yields  an  explicit  expression  for  4>i 
<&i(u7  vy  w)  = 

CwPd  (tiy  Vy  sgn (w)  \/w2+p2^ 

4n3/3vsIoSgn.(u>) \/w2+jP  Tf  ^sgn(w) y/ w2+p2^  R(-u7  -v) 

(17) 


At  last,  the  distribution  of  the  microwave  energy  deposition  can 
be  reconstructed  from  4>i  by  3-D  inverse  Fourier  transform. 
Equation  (17)  gives  an  exact  reconstruction  algorithm  for  planar 
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TAT  for  the  first  time.  Furthermore,  the  effects  of  the  finite  size 
of  the  detectors  and  the  finite  length  of  the  excitation  pulse  are 
included  explicitly.  From  (17),  it  can  be  inferred  that  the  re¬ 
constructed  image  spectrum  $d(u,  v,  w)  from  the  experimental 
data  without  the  consideration  of  these  two  effects,  as  was  pre¬ 
sented  by  previous  researchers  [9],  [20],  is  related  to  the  actual 
image  spectrum  3>i(u,  v,  w)  by 

4><*(w,  vy  w) 

=  47t2t/  ^ga(w) R(—u ,  — v,  w).  (18) 

Both  of  the  effects  result  in  multiplications  of  a  function  to  the 
actual  image  spectrum  in  the  frequency  domain.  They  are  equiv¬ 
alent  to  convolutions  in  the  spatial  domain,  which  blur  the  re¬ 
constructed  image.  However,  given  the  pulse  shape  and  the  sur¬ 
face  configuration  of  the  detector  surface,  the  two  effects  can  be 
reduced  by  deconvolution. 

To  summarize,  the  reconstruction  procedure  consists  of  the 
following  steps. 

1)  The  signal  from  the  detector  pd (x,  y ,  t)  is  Fourier  trans¬ 
formed  with  respect  to  t  to  yield  pd(x ,  y ,  k ).  Deconvolu¬ 
tion  with  respect  to  the  finite  pulse  length  can  be  imple¬ 
mented  immediately  after  the  Fourier  transform. 

2)  pd(: r,  t/,  k)  is  Fourier  transformed  with  respect  to  x  and 

2/,  yielding  Pd{u,  v,  *).  __ 

3)  According  to  (16)  and  (17),  Pd(u ,  v,  k)  is  mapped  to 

v ,  w). 

4)  Vy  «/)  is  deconvoluted  with  respect  to  the  finite  size 
of  the  detector,  giving  #i(u,  v,  w). 

5)  4>i  (u,  v ,  w)  is  inversely  Fourier  transformed  with  respect 
to  w,  v,  w  to  yield  <p(x',  ?/,  z'). 

The  order  of  steps  4)  and  5)  can  be  exchanged  so  that  more  stable 
deconvolution  algorithms  can  be  applied.  In  numerical  calcula¬ 
tions,  Pd(u ,  v,  k)  is  obtained  only  at  discrete  points;  hence  the 
mapping  from  Pd(u ,  v ,  k)  to  4>d(w,  v,  w)  needs  interpolation, 
which  can  be  a  major  source  of  distortion. 

B.  System  Setting 

The  experimental  setup  was  reported  in  [27]  and,  for  con¬ 
venience,  is  only  briefly  described  here  (Fig.  1).  The  x  axis 
points  perpendicularly  to  the  drawing  plane;  the  y  axis  points 
to  the  right  in  the  plane;  and  the  z  axis  points  downward  along 
the  acoustic  axis.  Microwave  pulses  are  transmitted  by  a  9-GHz 
microwave  generator.  The  pulse  width  is  0.5  ps.  The  object  to 
be  imaged  is  a  cylinder  of  pork  fat  containing  a  thin  layer  of 
connective  tissue  and  six  yellow  microstructures.  The  diameter 
of  the  cylinder  fat  is  14  mm  and  the  length  in  the  x  direction 
30  mm.  The  cylinder  was  immersed  in  mineral  oil  in  a  plexi¬ 
glass  tank.  The  central  frequency  of  the  ultrasonic  transducer 
(Panametrics)  is  2.25  MHz;  the  bandwidth  1 .8  MHz;  and  the 
diameter  of  the  active  element  6  mm.  More  details  about  the 
system  can  be  found  in  [27]. 

in.  Results  and  Discussion 

Our  method  was  applied  to  reconstructing  images  from  both 
the  simulated  and  the  experimental  data  in  a  two-dimensional 
(2-D)  case,  where  the  imaged  objects  were  uniform  along  the  x 
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Fig.  1.  Experimental  setup  for  TAT. 

axis.  Because  the  blur  due  to  the  finite  size  of  the  detector  sur¬ 
face  is  a  limiting  factor  on  the  resolution  of  images,  we  demon¬ 
strated  how  deconvolution  with  respect  to  the  detector  surface 
can  deblur  the  images.  We  chose  the  2-D  case  here  because  both 
the  computational  and  experimental  complexity  can  be  reduced 
more  in  the  2-D  case  than  in  the  3-D  one.  Nevertheless,  the  ex¬ 
tension  of  the  conclusions  of  the  2-D  case  to  the  3-D  one  is 
straightforward. 

A.  Simulation 

The  thermoacoustic  imaging  of  two  cylinders  was  numeri¬ 
cally  simulated.  Cylinders  were  chosen  because  the  analytical 
expression  for  their  thermoacoustic  signal  is  available  [28].  In 
the  simulations,  the  temporal-frequency  range  was  from  near 
0  to  1.5  MHz,  which  was  in  accordance  with  the  experimental 
one  and  with  our  previous  experiments  [11].  Two  simulations 
were  run.  The  first  one  was  to  test  our  reconstruction  algorithm 
under  an  ideal  experimental  condition,  which  is  noiseless  and 
does  not  consider  any  experimental  limitations  on  the  detectors. 
In  the  second  case,  the  effect  of  the  finite  size  of  the  detectors  on 
the  imaging  was  studied  while  noise  was  also  added.  Deconvo¬ 
lution  with  respect  to  the  finite  size  of  the  detector  surface  was 
applied  to  improve  the  lateral  resolution  of  the  blurred  image. 
Since  energy  deposition  is  a  positive  value,  only  the  positive 
components  of  the  reconstructed  image  were  retained,  and  the 
others  were  set  to  zero. 

In  step  3)  of  the  reconstruction,  which  is  the  mapping  from 
~Pd(u ,  v,  k)  to  $<*(«,  v ,  w ),  linear  interpolation  was  applied. 
By  adopting  the  zero-padding  technique  [25]  for  the  time-do¬ 
main  data,  one  can  increase  the  sampling  density  in  the  fc-space 
and,  consequently,  obtain  a  better  performance  of  the  interpo¬ 
lation  in  the  k- space.  In  the  reconstruction  from  the  simulation 
data  and  experimental  data,  we  appended  to  the  end  of  the  data 
the  same  number  of  zeros  as  in  the  original  collected  data,  so  that 
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Fig.  2.  The  reconstructed  image  of  the  cross-section  of  two  cylinders  with  a 
radius  of  2  mm  and  the  centers  separated  by  5.5  mm  under  ideal  experimental 
conditions. 

the  sampling  density  in  the  fc-space  was  doubled.  By  utilizing 
the  Wiener  filtering  method  [29],  deconvolution  with  respect  to 
the  finite  pulse  length  was  implemented  immediately  after  the 
Fourier  transform  with  respect  to  time  in  step  1).  As  the  decon¬ 
volution  with  respect  to  the  finite  size  of  the  detector  surface  is 
much  more  unstable  than  the  deconvolution  with  respect  to  the 
finite  pulse  length,  we  have  tried  two  methods  of  deconvolution: 
the  Wiener  filtering  method  and  the  piecewise  polynomial  trun¬ 
cated  singular  value  decomposition  (PP-TSVD)  [30]  method. 
The  first  method  can  be  implemented  in  the  spatial-frequency 
domain  and  is  more  computationally  efficient  than  the  second, 
but  the  performance  of  the  second  method  is  much  better,  as  it 
can  restore  sharp  boundaries  blurred  by  the  convolution  while 
avoiding  the  appearance  of  artificial  oscillations  in  an  unstable 
deconvolution.  Therefore,  we  adopted  the  PP-TSVD  method  to 
process  the  images.  Since  the  models  in  our  simulation  and  ex¬ 
periment  were  uniform  along  the  x  axis,  one-dimension  decon¬ 
volution  was  applied. 

Fig.  2  shows  the  reconstructed  image  from  the  simulated  data 
under  the  ideal  experimental  condition,  where  the  radius  of  the 
two  cylinders  was  2  mm;  the  distance  between  the  centers  of 
the  cylinders  was  5.5  mm;  the  centers  of  the  cylinders  were  po¬ 
sitioned  in  the  plane  of  z  =  10  mm;  the  scanning  range  of  the 
detector  along  the  y  axis  was  90  mm  with  a  step  size  of  0.5  mm; 
and  the  thermoacoustic  signals  were  sampled  for  40  /xs  at  a  sam¬ 
pling  rate  of  50  MHz.  The  reconstructed  image  is  in  good  agree¬ 
ment  with  the  real  objects,  whose  outlines  are  plotted  as  dotted 
circles  in  Fig.  2.  The  dimension  of  the  cylinders  is  3.75  mm 
along  the  2  direction  and  4.7  mm  along  the  y  direction.  The 
cylinder  is  a  little  deformed  laterally,  which  is  due  to  the  finite 
scanning  range  of  the  detector. 

Fig.  3  shows  the  images  before  and  after  deconvolution  with 
respect  to  the  finite  size  of  the  detector  surface  in  a  case  similar 
to  our  experimental  conditions.  The  noise  was  added  to  the  ther¬ 
moacoustic  signals,  and  the  signal-to-noise  ratio  (SNR)  was  50; 
the  diameter  of  the  detector  was  6  mm.  All  of  the  other  param¬ 
eters  were  the  same  as  those  in  the  first  case.  The  image  before 
deconvolution  is  shown  in  Fig.  3(a).  The  dimension  of  the  im¬ 
ages  of  the  objects  is  3.5  mm  along  the  y  axis,  which  agrees  well 
with  the  real  one,  4  mm.  However,  along  the  z  axis,  the  images 
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Fig.  3.  The  reconstructed  images  for  the  same  two  cylinders  as  in  Fig.  2  from 
noisy  data  (a)  before  and  (b)  after  the  deconvolution  with  respect  to  the  detector 
surface. 

of  the  two  cylinders  were  blurred  and  consequently  merged  into 
one,  which  is  predicted  by  our  analysis  of  the  effect  of  the  fi¬ 
nite  size  of  the  detector.  The  image  shows  no  clear  boundaries 
of  the  objects  along  the  y  axis.  After  deconvolution,  the  lateral 
boundaries  of  the  objects  become  very  clear  and  the  width  of 
the  objects  in  Fig.  3(b)  is  4. 1  mm,  which  is  quite  close  to  reality. 
Furthermore,  the  two  objects  can  be  distinguished  clearly.  After 
comparing  Fig.  3(a)  with  (b),  it  seems  that  the  ghost  images  be¬ 
come  slightly  more  obvious,  which  is  a  disadvantage  of  decon¬ 
volution.  Nevertheless,  it  is  obvious  that  deconvolution  with  re¬ 
spect  to  the  finite  size  of  the  detector  surface  can  improve  the 
lateral  resolution  greatly. 

In  Figs.  2  and  3,  there  are  some  ghost  images.  In  principle, 
our  reconstruction  method  is  exact  under  the  assumption  of 
thermal  confinement  and  constant  acoustic  speed.  However, 
several  factors  may  introduce  distortions.  First,  as_mentioned 
at  the  end  of  part  Section  II-A,  the  mapping  from  Pd(u ,  v,  k) 
to  $4(11,  v ,  w)  needs  interpolation,  which  is  a  major  source 
of  distortion.  This  distortion  can  be  reduced  by  increasing 
sampling  time  or  applying  a  better  interpolation  algorithm 
in  the  mapping.  Second,  in  experiments,  the  detector  cannot 
be  scanned  over  the  whole  plane.  Nevertheless,  Fig.  2  shows 
that  collecting  data  within  a  finite  area  of  the  collection  plane 
can  produce  images  of  sufficient  definition  to  determine  the 
configuration  and  position  of  the  objects. 

B.  Experimental  Result 

Fig.  4  shows  the  experimental  result.  The  images  before 
and  after  deconvolution  with  respect  to  the  finite  size  of  the 
detector  surface  are  shown  in  Fig.  4(a)  and  (b),  respectively. 
Fig.  4(c)  is  the  cross  section  of  the  biological  tissue,  which 
was  a  cylinder  with  a  radius  of  about  14  mm  and  3  cm  long.  It 
consisted  of  two  parts  of  fat  separated  by  a  very  thin  layer  of 
connective  tissue,  which  is  labeled  as  (7)  in  the  middle  of  the 
sample.  There  were  some  yellow  microstructures  among  the 
fat,  labeled  from  (l)-(6),  respectively.  Fig.  4(a)  is  the  image 
reconstructed  from  the  experimental  data  before  deconvolution. 
The  connective  tissue  between  the  two  parts  of  fat  and  the 
yellow  microstructures  are  imaged  clearly.  The  dimension 
of  the  image  is  16.4  mm  along  the  2  direction  and  19.2  mm 
along  the  y  direction.  However,  it  is  obvious  that  the  image 
before  deconvolution  is  blurred  along  the  y  axis,  which  makes 
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Fig.  4.  The  reconstructed  images  from  the  experimental  data  (a)  before  and 
(b)  after  the  deconvolution  with  respect  to  the  detector  surface;  (c)  the  cross 
section  of  a  cylinder  of  fat  sample  containing  six  yellow  microstructures  labeled 
from  (1H6)  and  a  layer  of  connective  tissue  in  the  middle  labeled  as  (7). 

the  lateral  boundaries  unclear  and  the  yellow  microstructures 
(1)  and  (2),  (3)  and  (4)  merge  into  one  object,  respectively. 
The  lateral  resolution  of  the  image  needs  to  be  improved. 
Consequently,  deconvolution  with  respect  to  the  finite  size  of 
the  detector  surface  was  applied  to  Fig.  4(a),  and  the  result  is 
shown  in  Fig.  4(b).  The  lateral  resolution  of  the  image  after 
deconvolution  is  much  improved.  The  merged  objects  can  be 
distinguished  clearly,  and  the  lateral  boundaries  of  the  cylinder 
become  much  clearer.  The  dimension  of  the  image  is  16.4  mm 
along  the  z  direction  and  16.7  mm  along  the  y  direction. 

C.  Discussion 

There  are  several  advantages  of  our  reconstruction  method. 
The  first  one  is  that  it  is  an  exact  reconstruction  algorithm. 
Unlike  other  reconstruction  methods  for  TAT  that  are  approxi¬ 
mate  ones,  our  reconstruction  method  provides  a  solid  base  for 
analyzing  and  improving  the  quality  of  reconstructed  images. 
Furthermore,  the  exact  reconstruction  method  has  a  broader 
application  than  the  approximate  ones.  For  example,  in  both 
our  simulation  and  experiment,  the  closest  distance  between  the 
objects  and  the  detectors  was  only  about  1  cm;  this  is  possible 
because  in  principle  there  is  no  limitation  on  the  detector-ob¬ 
ject  distance  in  our  method.  In  other  words,  the  detector  can 
be  placed  very  close  to  the  object  to  ensure  a  better  SNR. 
The  second  advantage  of  our  method  is  that  it  can  explicitly 
include  the  effect  of  many  limitations  from  the  experiment, 
such  as  the  finite  size  of  detector  surface,  the  microwave  pulse 


length,  and  the  finite  frequency  response  range  of  the  detector. 
Actually,  these  analyses  are  also  valid  for  other  approximate 
reconstruction  methods  as  long  as  the  other  reconstruction 
methods  are  able  to  produce  images  approximating  the  real 
objects.  Consequently,  our  analysis  of  the  blur  caused  by  the 
various  experimental  limitations  can  also  be  very  useful  for 
eliminating  the  limitations  in  other  reconstruction  methods. 
Lastly,  since  the  reconstruction  in  our  method  is  implemented 
in  the  frequency  domain,  the  efficiency  of  computation  is  much 
better  than  the  algorithm  implemented  in  the  time  domain  due 
to  the  use  of  the  efficient  Fourier  transformation  in  our  method. 
This  is  especially  important  for  real-time  3-D  imaging. 

From  the  above  images,  it  can  be  seen  that  there  is  no 
speckle  in  the  reconstructed  image.  Speckles  are  an  important 
factor  limiting  the  quality  of  pure  ultrasonic  imaging.  In  our 
technology,  the  detected  signals  are  directly  from  the  primary 
acoustic  waves  rather  than  reflective  or  scattered  waves. 
Further,  the  temporal  frequency  of  the  acoustic  signals  lies  in 
a  range  from  0  to  1.5  MHz,  which  is  only  weakly  scattered 
in  the  tissues.  The  above  two  factors  guarantee  that  there  is 
no  obvious  speckle  in  our  experimental  images.  However, 
the  issue  of  image  speckle  in  more  realistic  medical  imaging 
applications  of  our  algorithm  is  a  topic  for  future  consideration. 

The  formulas  in  this  paper  are  for  TAT  in  planar  geometry 
only.  However,  for  cylindrical  geometry  [19],  we  can  predict 
that  the  lateral  resolution  of  images  can  also  be  improved  by 
deconvolution  with  respect  to  the  detector  surface,  where  the 
deconvolution  is  carried  out  in  a  cylindrical  surface  instead  of 
a  plane.  For  spherical  geometry  [18],  similar  work  can  be  con¬ 
ducted  as  well. 

For  many  medical  imaging  applications,  the  acoustic  speed 
may  not  be  constant.  For  example,  the  acoustic  speed  inside  the 
female  breast  may  typically  exhibit  a  10%  variation;  however, 
our  simulation,  to  be  reported  elsewhere,  showed  that  the  image 
distortion  is  relatively  small. 

IV.  Conclusion 

We  developed  a  Fourier-domain  reconstruction  for  TAT  and 
obtained  an  exact  and  fast  reconstruction  algorithm.  The  effects 
of  the  finite  size  of  the  detectors  and  the  finite  length  of  the  ex¬ 
citation  pulse  were  included  explicitly  in  the  reconstruction  al¬ 
gorithm.  The  reconstruction  algorithm  was  verified  by  both  nu¬ 
merical  simulations  and  experimental  results.  Our  simulations 
demonstrated  that  the  blurring  caused  by  the  finite  size  of  the 
detector  surface,  which  is  a  key  limiting  factor  on  the  resolu¬ 
tion  of  images,  can  be  alleviated  by  deconvolution  with  respect 
to  the  detector  surface.  Other  effects  that  may  cause  the  blur  of 
the  images  can  be  treated  in  a  similar  way.  In  the  initial  exper¬ 
iment,  an  image  in  good  agreement  with  the  real  objects  was 
reconstructed  and  the  deconvolution  improved  the  resolution  of 
the  imaging  system.  The  method  can  also  be  extended  to  other 
configurations  of  data  collection. 
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Abstract — Microwave-induced  thermoacoustic  tomography 
(TAT)  in  a  cylindrical  configuration  Is  developed  to  image  bio¬ 
logical  tissue.  Thermoacoustic  signals  are  acquired  by  scanning  a 
flat  ultrasonic  transducer.  Using  a  new  expansion  of  a  spherical 
wave  in  cylindrical  coordinates,  we  apply  the  Fourier  and  Hankel 
transforms  to  TAT  and  obtain  an  exact  frequency-domain  recon¬ 
struction  method.  The  effect  of  discrete  spatial  sampling  on  image 
quality  is  analyzed.  An  aliasing-proof  reconstruction  method  is 
proposed.  Numerical  and  experimental  results  are  included. 

Index  Terms — Cylindrical,  frequency-domain  reconstruction, 
thermoacoustic  tomography. 


I.  Introduction 

THERMOACOUSTIC  tomography  (TAT)  combines  the 
strength  of  traditional  microwave  imaging  and  ultrasound 
imaging  [1]— [14].  Reviews  on  TAT  and  related  techniques 
can  be  found  in  [11],  [12],  [14].  Recently,  we  derived  exact 
reconstruction  algorithms  for  TAT  in  both  planar  and  spherical 
configurations;  these  are  reported  in  the  companion  papers  [1 1], 
[12].  We  recognize,  however,  that  in  some  applications  such  as 
the  imaging  of  the  limbs,  a  cylindrical  scanning  surface  may  be 
more  appropriate.  In  this  paper,  using  a  new  expansion  formula 
in  cylindrical  coordinates,  we  derive  a  frequency-domain 
reconstruction  algorithm  [15]-[19]  and  report  our  numerical 
and  experimental  results  in  two-dimensional  (2-D)  cases. 

n.  Methods 

We  assume  that  the  detector  scans  on  a  cylindrical  surface 
with  a  radius  of  p,  which  encircles  all  microwave  absorbing  ob¬ 
jects.  In  our  paper,  a  coordinate  with  a  prime  refers  to  the  po¬ 
sition  in  an  imaged  object,  while  a  coordinate  without  a  prime 
refers  to  that  of  a  detector.  In  the  case  of  thermal  confinement. 
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the  temporal  spectrum  of  acoustic  field  p(r,  k)  is  related  to  the 
microwave  absorption  distribution  <p( r')  by  the  following  equa¬ 
tion  [11]: 

Hr,  k)  -  JjJ  gHfc:  ^1)  dr’ 

where  the  symbols  are  defined  as  in  [11].  Cylindrical  coordi¬ 
nates  are  used  in  the  following  derivation,  where  £  is  shown 
in  [12,  Fig.  2],  and  p,  <j>  are  the  polar  coordinates  within  the 
x-y  plane.  Following  the  derivation  of  the  series  expansion  of 
l/|r  -  r'|  [20],  we  obtained  the  following  new  identity  for  a 
series  expansion  of  a  spherical  wave  in  a  cylindrical  coordinate 
system  (see  the  Appendix  for  the  derivation): 


exp(— ik\r  —  r'[) 

47r|r  —  r'| 

-i 

=  —  /  dkz  exp[— ikz(z*  —  z)\ 
7-00 


oo 

•  A(m,  pp\  pp)  exp -  </>)] 

m— — oo 


(2) 


where  p  =  sgn(fc)  ^ |/ c2  -  fcff;  sgn( )  is  the  signum  function; 
and  A  is  the  function  defined  as 

( JM)Hl(pp),  if|*|>|M 
A(m,  pp\  /*)  =  <»  JMf/)Krn(\p\p),  if  \k\  <  \kz\ 

\  7 r 

where  Jmt  H^,  Jm,  and  Km  are  the  rath-order  Bessel, 
second-kind  Hankel,  and  modified  Bessel  functions,  respec¬ 
tively.  It  has  been  assumed  in  the  above  two  equations  that 
p  >  p'.  Substituting  (2)  into  (1)  results  in 

/OO 

dkz  exp [—ikz{zt  —  z)\ 

-oo 

oo 

•  A(m,  pp' ,  pp)  exp[-im(cj)' -  <f>)].  (3) 

m=— oo 

The  | k |  >  \kz\  part  of  the  integration  with  respect  to  kz  rep¬ 
resents  the  contribution  from  the  propagation  wave,  while  the 
\k\  <  \kz\  part  represents  the  evanescent  wave.  As  the  evanes¬ 
cent  wave  decays  rapidly  at  a  distance  several  wavelengths  from 
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the  source,  it  is  not  suitable  for  thermoacoustic  imaging.  For  the 
case  of  |Jfc|  >  1^1,  after  Fourier  transforming  both  sides  of  the 
above  equation  with  respect  to  (j>  and  z,  we  have 


Pi(m,  kz,k)  = 


8tt  C 


dp'p'pi(m,  kz,  p')Jm(pp') 


(4) 


where  Pi(m,  kz,  k)  and  kz,  p')  are  the  Fourier  trans¬ 

forms  of  p(r,  k)  and  ip( r'),  respectively.  Noticing  that  the  right 
side  of  (4)  is  actually  a  Hankel  transform,  an  inverse  Hankel 
transform  gives 


<pi(m,  k,, 


r 

’  PvJo  Jo 


PPi  (m,  k^,  k)Jln{pp') 

kfj(k)H?M  ’ 
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Applying  a  variable  change  of  the  integral  variable  from  p  to  k 
to  the  above  equation  results  in 


<Pi{m,  k. 
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pi (777,  kz:  ) 

;  rj(k)HMw)  ’ 


|*|  >  | kz 


(5) 


At  last,  kz ,  p')  is  inversely  Fourier  transformed  with  re¬ 

spect  to  m  and  kz  to  yield  <p(</>',  z\  p').  Equation  (5)  gives  an 
exact  mapping  relation  between  the  spectrum  of  the  collected 
signals  and  the  spectrum  of  the  distribution  of  microwave  en¬ 
ergy  deposition  and  is  the  essence  of  our  reconstruction  method. 

An  exact  reconstruction  method  for  ultrasonic  reflectivity 
imaging  with  a  cylindrical  scanning  surface  was  given  in  [16]. 
However,  our  results  are  much  simpler  and  more  stable.  In  their 
equation  A24,  7m(/rr0),  where  r0  is  the  radius  of  the  scanning 
cylindrical  surface,  appeared  in  the  denominator  and  can  be 
zero  for  some  values  of  p;  consequently,  this  term  can  cause 
instability.  In  our  (5),  H^(pp)  appeared  in  the  denominator, 
which  cannot  be  zero  for  a  finite  p. 

To  summarize,  the  reconstruction  procedure  consists  of  the 
following  steps. 

1)  The  signal  from  the  detector  p(</>,  z,  t)  is  Fourier  trans¬ 
formed  with  respect  to  t  to  yield  p(<j> ,  z,  k).  Deconvolu¬ 
tion  with  respect  to  the  finite  pulse  length  can  be  imple¬ 
mented  immediately  after  the  Fourier  transform. 

2)  p(</>,  z,  k)  is  Fourier  transformed  with  respect  to  *  and  0, 
giving  Pi  (m,  kxy  k). 

3)  According  to  (5),  px(m,  kz,  k)  is  mapped  to 
<pi(m,  kZJ  p'). 

4)  (pi  (m,  kz,  p')  is  inversely  Fourier  transformed  with  re¬ 
spect  to  m,  kz  to  yield  z! ,  pf). 


III.  Results  and  Discussion 

To  test  our  method,  images  from  both  numerically  simulated 
and  experimental  data  were  reconstructed  in  a  2-D  case.  We 
chose  the  2-D  case  rather  than  the  three-dimensional  (3-D)  case 
to  reduce  the  computational  and  experimental  complexity.  For 
the  2-D  case,  the  reconstruction  equation  can  be  derived  from 


(a)  (b) 
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Fig.  1.  The  images  (a)  before  and  (b)  after  the  reconstruction  from  the 
simulated  data  of  two  cylinders. 


Fig.  2.  (a)  The  cross  section  of  a  fat  sample  containing  5  pieces  of  muscle 

cylinders,  (b)  The  reconstructed  image  from  the  experimental  data, 

(4)  by  replacing  all  kz  with  zero.  The  extension  of  the  conclu¬ 
sions  of  the  2-D  case  to  a  3-D  one  is  straightforward. 

A.  Numerical  Simulation 

The  thermoacoustic  imaging  of  two  cylinders  was  numeri¬ 
cally  simulated,  where  the  radius  of  each  cylinder  was  2  mm; 
the  distance  between  the  centers  of  the  cylinders  was  5  mm;  and 
the  center  of  one  of  the  cylinders  was  positioned  at  the  origin  of 
the  circle  of  detection.  Cylinders  were  chosen  because  the  an¬ 
alytical  expression  for  their  thermoacoustic  signal  is  available 
[21].  In  the  simulations,  the  temporal-frequency  range  was  from 
about  0  to  2  MHz,  which  was  close  to  our  experimental  situa¬ 
tion  [14].  For  the  noiseless  simulated  data,  the  reconstruction  is 
almost  perfect.  Therefore,  we  show  only  the  results  from  noisy 
data.  Fig.  1  shows  the  images  before  and  after  the  reconstruc¬ 
tion  from  the  simulated  data  with  introduced  additive  noise.  The 
units  for  the  signals  and  energy  deposition  in  Figs.  1  and  2  are 
relative  ones.  Calibration  of  our  system  is  needed  to  obtain  an 
absolute  measurement.  The  radius  of  the  circle  of  detection  was 
30  mm;  the  angular  scanning  range  was  2tx  with  256  steps;  and 
the  thermoacoustic  signals  were  sampled  for  50  ps  at  a  sam¬ 
pling  rate  of  4  MHz.  The  signal-noise-ratio  (SNR)  of  the  raw 
data  shown  in  Fig.  1(a)  was  1.  The  reconstructed  image  shown 
in  Fig.  1(b)  is  in  good  agreement  with  the  real  objects,  whose 
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outlines  are  plotted  as  dotted  circles  in  Fig.  1(b).  The  dimen¬ 
sions  of  the  reconstructed  cylinders  are  4  mm  along  both  the 
x  and  the  y  directions.  The  SNR  of  the  reconstructed  image  is 
about  8,  which  is  improved  greatly  compared  with  that  of  the 
raw  data. 

B.  Experiment  Results 

The  experimental  setup  for  2-D  TAT  in  a  cylindrical  con¬ 
figuration  is  the  same  as  that  in  [12].  The  sample  is  shown  in 
Fig.  2(a),  which  was  photographed  after  the  experiment.  Mi¬ 
crowave  pulses  were  delivered  to  the  sample  from  below.  The 
imaging  plane  was  2  cm  above  the  bottom  of  the  tissue  sample. 
Above  the  plane,  there  is  another  layer  of  fat  about  1  cm  thick. 
The  sample  consisted  of  five  muscle  cylinders  with  a  diameter 
of  about  3  mm  and  a  height  of  6  mm.  The  muscle  cylinders  were 
surrounded  by  pork  fat.  The  electrical  property  of  interest  to  this 
imaging  technique  is  the  microwave  attenuation  coefficient  of 
the  medium  at  the  experimental  microwave  frequency,  3  GHz. 
The  microwave  attenuation  coefficients  of  fat  and  muscle  are 
9  cm”1  and  1  cm”1,  respectively.  The  microwave  absorption  in 
mineral  oil  can  be  neglected,  compared  with  the  absorption  in 
fat  and  muscle.  During  the  experiment,  the  transducer  scanned 
around  the  sample  at  a  radius  of  7.1  cm  from  0°  to  360°  with  a 
step  size  of  2.25°.  The  thermoacoustic  signals  were  sampled  for 
60  ps  at  a  sampling  rate  of  20  MHz.  The  time  between  the  end 
of  a  microwave  pulse  and  the  acquisition  of  the  thermoacoustic 
signal  was  between  10  ps  and  20  ps  in  our  system,  depending 
on  the  distance  of  the  transducer  to  the  nearest  sample  surface. 

Fig.  2(b)  shows  the  reconstructed  image  from  the  experi¬ 
mental  data.  The  reconstructed  image  is  in  good  agreement 
with  the  real  objects.  The  boundaries  between  the  fat  and  the 
surrounding  medium  and  the  muscle  cylinders  are  imaged 
clearly.  However,  it  can  be  seen  that  the  quality  of  the  image 
decreases  with  the  increasing  distance  of  the  objects  from  the 
center  of  the  circle  of  detection.  One  possible  reason  is  that  the 
finite  surface  area  of  the  detector,  which  has  a  6-mm  diameter  in 
this  experiment,  may  cause  blurring  of  the  image  perpendicular 
to  the  radial  direction,  and  this  blurring  is  more  serious  when 
the  object  is  farther  from  the  center.  Another  possible  reason  is 
that  the  microwave  field  decreases  when  the  radius  increases 
in  our  irradiation  configuration. 

Our  method  can  be  applied  to  analyze  the  effect  of  the  dis¬ 
crete  sampling  by  the  detector  along  the  circle  of  detection  on 
imaging.  This  can  be  illustrated  by  analyzing  the  signals  from  a 
point  source  located  at  radius  p\.  According  to  (4) 

Pi (wi,  k)  oc  Jm(fcpi).  (6) 

Fig.  3  shows  how  Jm{kpi)  changes  with  m,  where  k  — 
8.37  mm”1  (the  wave  number  of  a  2-MHz  acoustic  wave)  and 
Pi  =  10  mm.  It  is  clear  that  Jm(kpi)  has  considerable  value 
until  m  «  kpu  where  the  Bessel  function  makes  a  transition 
from  near-field  behavior  to  far-field  behavior.  Therefore,  it 
is  safe  to  claim  that,  with  respect  to  variable  <f>,  p( r,  k)  is 
band-limited  by  kpi.  According  to  the  Nyquist  criteria,  the 
number  of  scanning  points  per  cycle  should  be  at  least  2kpi  to 
avoid  aliasing.  In  other  words,  for  a  fixed  number  of  scanning 
points  JV,  the  maximum  wave  number  before  aliasing  occurs 


Fig.  3.  im(fcpi)  versus  m,  where  k  =  8.37  mm  1  (the  wave  number  of  a 
2  MHz  acoustic  wave)  and  pi  =  10  mm. 

is  fcmax  «  N/(2pi).  It  can  be  seen  that  the  maximum  wave 
number  is  inversely  proportional  to  pi.  For  the  same  N  and 
temporal  spectrum  of  signal,  the  aliasing  may  be  more  serious 
for  signals  coming  from  sources  at  a  greater  radial  distance 
than  for  those  closer  to  the  center.  The  above  analysis  also 
points  out  a  way  to  produce  an  aliasing-free  image  from  the 
data  obtained  by  discrete  detection.  That  is  to  apply  a  filter  in 
the  temporal-frequency  domain  to  the  spectrum  of  the  temporal 
data  with  a  stopband  at  about  JV/(2pmax),  where  pmax  is  the 
maximum  radius  of  imaging  range  of  interest.  The  application 
of  the  filter  will  decrease  the  resolution  of  the  image;  however, 
it  can  guarantee  that  there  will  be  no  aliasing  in  the  image. 

C.  Discussion 

Since  our  method  is  implemented  in  the  frequency  domain 
using  the  fast  Fourier  transform  (FFT)  technique,  the  computa¬ 
tional  efficiency  is  much  greater  than  if  implemented  in  the  time 
domain.  The  most  time-consuming  computation  in  the  numer¬ 
ical  reconstruction  lies  in  (5),  which  is  a  Hankel  transform.  For¬ 
tunately,  a  quasi-fast  algorithm  for  it,  which  is  as  efficient  as  a 
one-dimensional  FFT,  is  available  [22].  Following  the  methods 
in  [1 1],  our  method  can  explicitly  include  and  further  eliminate 
the  effect  of  many  limitations  from  the  experiment,  such  as  the 
finite  size  of  the  detector  surface,  the  microwave  pulse  length, 
and  the  finite  response  frequency  range  of  the  detector.  Addi¬ 
tionally,  combining  our  method  and  the  techniques  in  [16],  a 
new  exact  reconstruction  algorithm  for  3-D  ultrasonic  reflec¬ 
tivity  imaging  with  a  cylindrical  aperture  can  be  derived.  Finally, 
we  would  like  to  point  out  that  the  reconstruction  methods  re¬ 
ported  in  this  paper  and  the  two  companion  papers  [  1 1  ],  [  1 2]  are 
also  applicable  to  photoacoustic  or  optoacoustic  tomography  as 
well  as  other  diffraction-based  inverse  source  problem. 

The  size  of  tissue  samples  that  can  be  imaged  by  our  system 
is  mainly  limited  by  the  safety  standard  on  microwave  power, 
the  microwave  frequency,  the  microwave  irradiation  configu¬ 
ration,  the  sensitivity  of  the  ultrasonic  transducer,  the  dynamic 
range  of  the  preamplifier  and  sampling  system,  and  the  afford¬ 
able  imaging  time.  The  effect  of  microwave  frequency  on  the 
imaging  depth  was  addressed  in  reference  [13].  A  microwave 
irradiation  configuration  that  renders  a  uniform  microwave  irra¬ 
diation  within  the  sample  will  also  increase  the  capacity  of  the 
system  to  image  larger  samples.  A  large  dynamic  range  of  the 


832 


IEEE  TRANSACTIONS  ON  MEDICAL  IMAGING,  VOL.  21,  NO.  7,  JULY  2002 


preamplifier  and  the  sampling  system  is  necessary  to  accurately 
collect  the  thermoacoustic  signals  from  both  the  surface  and  the 
inside  of  a  sample.  A  more  sensitive  ultrasonic  transducer  and 
a  longer  imaging  time  can  improve  the  signal-to-noise  ratio  of 
acoustic  signals  and  make  the  weak  signals  from  the  inside  of 
large  samples  detectable. 

In  our  initial  computation,  the  reconstruction  of  a  single  2D 
image  required  about  2  min  in  a  Dell  Precision  330  computer 
(Intel  Pentium  4  processor  with  a  clock  frequency  of  1 .5  GHz) 
with  Matlab  programs  if  there  was  no  precomputation  of  Bessel 
and  Hankel  functions.  However,  our  initial  computation  was 
aimed  at  verifying  the  proposed  algorithm  rather  than  demon¬ 
strating  the  computation  efficiency.  The  proposed  algorithm  can 
be  implemented  with  high  computational  efficiency  as  stated  in 
the  discussion  section.  For  high  computational  efficiency,  the 
program  should  be  coded  with  languages  such  as  C  or  Fortran, 
Bessel  and  Hankel  functions  should  be  precomputed,  and  the 
fast  Hankel  transform  algorithm  should  be  adopted.  The  evalu¬ 
ation  of  the  computation  efficiency  of  our  algorithm  is  a  topic 
for  future  studies. 


IV.  Conclusion 

Using  a  new  expansion  of  a  spherical  wave  in  the  cylindrical 
coordinate  system,  we  applied  the  Fourier  transform  and  Hankel 
transform  techniques  to  TAT  with  a  cylindrical  detection  sur¬ 
face.  The  reconstruction  algorithm  is  verified  by  both  numerical 
simulations  and  experimental  results  in  2-D  cases.  The  method 
was  applied  to  analyze  the  effect  of  discrete  sampling  by  the  de¬ 
tector  along  the  circle  of  detection  on  imaging;  an  aliasing-free 
reconstruction  method  for  discrete  sampling  along  the  azimuth 
direction  is  proposed. 


Appendix 

The  derivation  of  (2)  will  be  presented  here.  The  spherical 
wave  Gk(ry  r')  =  exp(— iA;|r  —  r'|)/(47r|r  —  r'|)  is  a  solution 
to  the  wave  equation  with  a  point  source 

V^(?fc(r,  r')  +  k2Gh( r,  r')  =  -6( r  -  r').  (Al) 

The  solution  can  be  expanded  in  terms  of  orthonormal  functions 
of  z  and  (j>  in  a  cylindrical  coordinate  system 

1  \  2  oo  oo 

2  )  ^ ^  I  dkzgTH(ky  kzy  py  p  ) 

/  m=— oo  J  —  00 

•  exp[im(<f>  -  <j>)  *f  ikz(z  -  z')].  (A2) 


Substituting  (A2)  into  (Al)  results  in  an  equation  for  the  radial 
Green’s  function  gm 


1  d_ 
P  dp 


8{p  -  p') 

gm  = - — * 

(A3) 


When  \k\  <  |fc,|,  following  the  derivation  of  the  series  expan¬ 
sion  of  l/|r  -  r'|  [20],  one  obtains  a  similar  expansion  for  the 
spherical  wave 


(A4) 


We  next  consider  the  case  of  |  Ar|  >  \kz\  and  k  >  0.  Noticing  that 
when  p  — *  oo,  gni  behaves  asymptotically  as  exp [-ip(p  -  p')\ 
(p  >  p '  is  implicit  in  our  model),  one  can  follow  the  derivation 
in  [20]  and  obtain 

=  (A5) 

Similarly,  for  \k\  >  \kz\  and  k  <  0 

gm  =  ^Jm{WWMp)-  (A6) 

Using  the  following  identities  of  Bessel  and  Hankel  functions 
[23]: 

Jm{PP)=(~  1)  pp) 

and  combining  (A2)  and  (A4)-(A6),  we  obtain  (2). 
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Abstract  A  study  of  pulsed-microwave-induced  thermoacoustic  tomography  in  biological  tissues  is 
presented.  A  backprojection  algorithm  based  on  rigorous  theory  is  used  to  reconstruct  the  cross-sectional 
image  from  the  thermoacoustic  measurement  in  a  circular  configuration  that  encloses  the  sample  under 
study.  The  results  demonstrate  the  possibility  of  application  in  detecting  small  tumors  buried  in  biological 
tissues  using  microwave  absorption  contrast  and  ultrasound  spatial  resolution.  Finally,  the  method  is 
compared  with  laser-induced  thermoacoustic  tomography. 
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OCIS  Codes:  (170.5120)  Photoacoustic  imaging;  (170.6960)  Tomography 

1.  Introduction 

Pulsed  microwave-induced  thermoacoustic  tomography  combines  the  advantages  of  both  ultrasound  spatial 
resolution  and  microwave  absorption  contrast  [l]-[4].  With  this  technique,  a  very  short  microwave  pulse  (<1 
microsecond)  heats  a  sample;  the  sample  then  absorbs  the  microwave  energy  in  a  confined  time  and  simultaneously 
generates  temporal  thermoacoustic  waves,  which  are  strongly  related  to  the  locally  absorbed  microwave  energy.  The 
thermoacoustic  signals  have  a  wide  frequency  range  up  to  -  MHz  and  carry  the  information  about  microwave 
absorption  distribution  with  millimeter  spatial  resolution.  In  practice,  microwaves  at  300  MHz  ~  3  GHz  with  0.1  ~  1 
(is  pulse  are  often  adopted,  which  provide  several  centimeters  penetration  depths  in  biological  tissues.  Due  to  the 
bounded  water  and  salt  that  exist  in  cancer  cells,  a  tumor  absorbs  more  microwave  energy  and  generates  more 
intense  thermoacoustic  waves  than  the  surrounding  tissues  [5],  [6].  The  wide  range  of  absorption  values  among 
various  tissues  makes  it  possible  to  achieve  a  high  image  contrast.  In  addition,  the  long  penetration  depth  allows  this 
technique  to  detect  interior  tumors. 

In  this  paper,  we  present  our  study  of  pulsed-microwave-induced  thermoacoustic  tomography  under  a  circular 
measurement  configuration  in  biological  tissues.  A  wide  beam  of  short-pulse  microwave  energy  is  used  to  illuminate 
a  sample  from  the  bottom.  An  unfocused  ultrasonic  transducer  with  a  small  aperture  is  used  to  record  the 
thermoacoustic  signals.  A  backprojection  method  based  on  rigorous  theory  is  used  to  reconstruct  the  cross-sectional 
image  from  the  measured  data.  A  phantom  sample  is  investigated.  The  reconstructed  image  agrees  with  the  original 
sample  very  well.  Finally,  the  method  is  compared  with  laser-induced  thermoacoustic  tomography. 

2.  Method 


Fig.  1  Schematic  of  the  circular  measurement  system 


A  schematic  view  of  the  circular  measurement 
system  for  our  study  is  shown  in  Fig.  1.  A 
plexiglass  container  is  filled  with  mineral  oil. 
An  unfocused  transducer  is  immersed  inside  it 
and  fixed  on  a  rotation  device.  A  step  motor 
drives  the  rotation  device  and  then  moves  the 
transducer  scan  around  the  sample  on  a 
horizontal  x-y  plane,  where  the  transducer 
horizontally  points  to  the  rotation  center.  A 
sample  is  immersed  inside  the  container  and 
placed  on  a  holder,  which  is  made  of  thin 
plastic  material  that  is  transparent  to  the 
microwave.  The  transducer  (V323, 
Panametrics)  has  a  central  frequency  of  2.25 
MHz,  and  a  diameter  of  6mm.  The  microwave 
pulses  transmitted  from  a  3-GHz  microwave 
generator  have  a  pulse  energy  of  10  mJ  and  a 
pulse  width  of  0.5ps.  A  function  generator 
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(Protek,  B-180)  is  used  to  trigger  the  microwave  generator,  control  its  pulse  repetition  frequency,  and  synchronize 
the  oscilloscope  sampling.  Microwave  energy  is  delivered  to  the  sample  by  a  rectangular  waveguide  with  a  cross 
section  of  72  mm  x  34  mm.  A  personal  computer  is  used  to  control  the  step.  The  signal  from  the  transducer  is  first 
amplified  through  a  pulse  amplifier,  then  recorded  and  averaged  500  times  by  an  oscilloscope  (TDS640A, 
Tektronix),  and  finally  transferred  to  a  personal  computer  for  imaging. 

We  assume  the  tissue  to  have  inhomogeneous  microwave  absorption  and  a  relatively  homogeneous  acoustic 
property.  When  the  microwave  pulse  duration  is  <  1  (is,  the  heat  diffusion’s  effect  on  the  thermoacoustic  wave  in  the 
tissue  can  be  ignored.  The  speed  of  sound  c  in  most  soft  tissue  is  relatively  constant  at  ~1.5  mm/ps.  Suppose  a  delta 
illuminating  function  /0§(/)  and  a  detected  acoustic  pressure  p(r0,t)  on  the  circular  surface  r  =  r0  =(p0,(p0,z0)  and 

time  t.  If  challenged  to  detect  small  size  tumors,  we  can  safely  remove  the  low-frequency  component.  In  addition, 
the  wavelengths  of  the  high-frequency  thermoacoustic  waves  are  much  smaller  than  the  detecting  distances  between 
the  thermoacoustic  sources  and  the  transducers.  Under  the  above  conditions,  the  spatial  absorption  function  A( r) 
can  be  calculated  by  the  following  2D  surface  integral  in  the  cylindrical  configuration 


where  p  is  the  isobaric  volume  expansion  coefficient  and  Cp  is  the  heat  capacity.  The  above  reconstruction  formula 


(z0-z)2 1  dp(r09t) 
I  r-r 0\2  t  dt 


indicates  that  the  cross-sectional  image  of  any  z  plane  is  determined  mainly  by  the  data  measured  on  the  circle  of  the 
same  z  plane.  In  other  words,  if  small  absorption  sources  are  located  on  a  z  plane,  a  set  of  circular  measurement  data 
on  the  same  plane  can  be  sufficient  to  yield  a  good  cross-sectional  image. 


3.  Results  and  Discussion 


A  phantom  sample  was  imaged  by  our  microwave-induced  thermoacoustic  tomography  system.  The  measurement 
diagram  and  its  cross-sectional  photograph  are  shown  in  Fig.  2  (a)  and  (b),  respectively.  Three  small  absorbers, 
which  were  made  of  gelatin,  salt  and  water,  were  buried  inside  a  large  fat  base.  The  transducer  rotationally  scanned 
the  sample  from  0  to  360  degrees  with  a  step  size  of  2.25  degree.  The  reconstructed  image  produced  by  our 
backprojection  method  is  shown  in  Fig.  2  (c),  which  agrees  with  the  original  sample  very  well.  The  relative  locations 
and  sizes  of  those  thermoacoustic  sources  perfectly  match  the  buried  objects  in  the  original  sample. 


Fig.  2  (a)  Diagram  of  the  measurement  scheme,  (b)  Photograph  of  the  sample,  (c)  Reconstructed  image. 


Lihong  Wang  —  Appendix  4  (Page  3  of  3) 


x  (mm) 


(a) 


0.04  mm 

— HK— 


0.19  mm 


0.23  mm 


(b) 


E 

E 


0.0  0.2  0.4  0.6  0.8  1.0 


Fig.  3  (a)  The  schematic  of  two  adjacent  dots,  (b)  The  reconstructed  image. 

For  comparison,  as  Fig.  3  shows,  an  example  of  two  black  dots  with  a  gap  distance  40  |im  was  imaged  by  traditional 
photoacoustic  tomography,  in  which  short  laser  pulses  were  used  to  illuminate  the  sample.  A  10  MHz  transducer 
was  used  to  detect  photoacoustic  signal.  Because  the  laser  duration  is  only  4.7  ns,  a  higher  frequency  thermoacoustic 
signal  was  generated  and  a  higher  spatial  resolution  was  obtained.  However,  the  laser  was  not  able  to  penetrate  very 
deeply  inside  of  the  tissue. 

3.  Conclusion 

We  have  presented  our  study  on  pul sed-micro wave-induced  thermoacoustic  tomography  by  a  circular  measurement 
configuration  in  biological  tissues.  A  backprojection  algorithm  is  used  to  reconstruct  the  cross-sectional  images.  The 
reconstructed  image  of  a  phantom  sample  agrees  with  the  original  values  very  well.  The  result  demonstrates  that  the 
circular  measurement  configuration  combined  with  the  backprojection  method  is  a  promising  technique  for  use  in 
detecting  small  tumors  buried  in  biological  tissues  with  microwave  absorption  contrast  and  ultrasound  spatial 
resolution  (~  mm).  Alternatively,  a  pulsed-laser  technique,  which  provides  laser  absorption  contrast  as  well  as  high 
spatial  resolution  (~  50  Jim),  can  be  employed  as  an  illumination  source  to  detect  small  tumors  that  are  not  very 
deeply  buried. 

Acknowledgments 

This  project  was  sponsored  in  part  by  the  U.S.  Army  Medical  Research  and  Materiel  Command  Grant  No. 
DAMD 17-00- 1-0455,  the  National  Institutes  of  Health  Grant  No.  R01  CA71980,  the  National  Science  Foudation 
Grant  No.  BES-9734491,  and  Texas  Higher  Education  Coordinating  Board  Grant  No.  ARP  000512-0123-1999. 

References 

1  R.  A.  Kruger,  K.  K.  Kopecky,  A.  M.  Aisen,  D.  R.  Reinecke,  G.  A.  Kruger  and  W.  L.  Kiser,  "Thermoacoustic  CT  with  radio 
waves:  A  medical  imaging  paradigm,"  Radiology  211,  pp.  275-278,  1999. 

2  G.  Ku  and  L.-H.  V.  Wang,  "Scanning  thermoacoustic  tomography  in  biological  tissues,"  Med.  Phys.  27,  pp.  1195-1202, 
2000. 

3  G.  Ku  and  L.-H.  V.  Wang,  "Scanning  microwave-induced  thermoacoustic  tomography:  Signal,  resolution,  and  contrast," 
Med.  Phys.  28,  pp.  4-10, 2001. 

4  M.  H.  Xu  and  G.  Ku,  and  L.-H.  V.  Wang,  "Microwave-induced  thermoacoustic  tomography  using  multi-sector  scanning", 
Med.  Phys.  29,  pp.  1958-1963, 2001. 

5  S.  Chaudhary,  R.  Mishra,  A.  Swarup,  and  J.  Thomas,  "Dielectric  properties  of  normal  human  breast  tissues  at  radiowave  and 
microwave  frequencies,"  Indian  Journal  of  Biochemistry  and  Biophysics  21, 76-79  (1984). 

6  W.  Joines,  Y.  Zhang,  C.  Li,  and  R.  Jirtle,  "The  measured  electrical  properties  of  normal  and  malignant  human  tissues  from 
50-900  MHz,"  Medical  Physics  21,  547-550  (1994). 


